{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cathedral-trance",
   "metadata": {},
   "source": [
    "# Jelly roll model\n",
    "\n",
    "In this notebook we show how to set up and solve the \"two-potential\" model from \"Homogenisation of spirally-wound high-contrast layered materials\", S. Psaltis, R. Timms, C.P. Please, S.J. Chapman, SIAM Journal on Applied Mathematics, 2020.\n",
    "\n",
    "We consider a spirally-wound cell, such as the common 18650 lithium-ion cell. In practice these cells are constructed by rolling a sandwich of layers containing the active cathode, positive current collector, active cathode, separator, active anode, negative current collector, active anode, and separator. The \"two-potential\" model consists of an equation for the potential $\\phi^\\pm$ in each current collector. The potential difference drives a current $I$ through the electrode/separator/electrode sandwich (which we refer to as the \"active material\" in the original paper). Thus, in non-dimensional form, the model is \n",
    "\n",
    "$$ \\frac{\\delta^+\\sigma^+}{2\\pi^2}\\frac{1}{r}\\frac{\\mathrm{d}}{\\mathrm{d}r}\\left(\\frac{1}{r}\\frac{\\mathrm{d}\\phi^+}{\\mathrm{d}r}\\right) + 2I(\\phi^+-\\phi^-) = 0,$$\n",
    "$$ \\frac{\\delta^-\\sigma^-}{2\\pi^2}\\frac{1}{r}\\frac{\\mathrm{d}}{\\mathrm{d}r}\\left(\\frac{1}{r}\\frac{\\mathrm{d}\\phi^-}{\\mathrm{d}r}\\right) - 2I(\\phi^+-\\phi^-) = 0,$$\n",
    "with boundary conditions \n",
    "$$ \\frac{\\mathrm{d}\\phi^+}{\\mathrm{d}r}(r=r_0) = 0, \\quad \\phi^+(r=1) = 1, \\quad \\phi^-(r=0) = 0, \\quad \\frac{\\mathrm{d}\\phi^-}{\\mathrm{d}r}(r=1) = 0.$$\n",
    "\n",
    "For a complete description of the model and parameters, please refer to the original paper.\n",
    "\n",
    "It can be shown that the active material can be modelled using any 1D battery model we like to describe the electrochemical/thermal behaviour in the electrode/separator/electrode sandwich. Such functionality will be added to PyBaMM in a future release and will enable efficient simulations of jelly roll cells.  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "whole-diabetes",
   "metadata": {},
   "source": [
    "## Two-potential resistor model\n",
    "In this section we consider a simplified model in which we ignore the details of the anode, cathode and separator, and treat them as a single region of active material, modelled as an Ohmic conductor, with two such regions per winding. In this case the model becomes \n",
    "\n",
    "$$ \\frac{\\delta^+\\sigma^+}{2\\pi^2}\\frac{1}{r}\\frac{\\mathrm{d}}{\\mathrm{d}r}\\left(\\frac{1}{r}\\frac{\\mathrm{d}\\phi^+}{\\mathrm{d}r}\\right) + \\frac{2\\sigma^{a}(\\phi^--\\phi^+)}{l\\epsilon^4} = 0,$$\n",
    "$$ \\frac{\\delta^-\\sigma^-}{2\\pi^2}\\frac{1}{r}\\frac{\\mathrm{d}}{\\mathrm{d}r}\\left(\\frac{1}{r}\\frac{\\mathrm{d}\\phi^-}{\\mathrm{d}r}\\right) - \\frac{2\\sigma^{a}(\\phi^--\\phi^+)}{l\\epsilon^4} = 0,$$\n",
    "along with the same boundary conditions.\n",
    "\n",
    "We begin by importing PyBaMM along with some other useful packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "heard-cartridge",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m A new release of pip is available: \u001b[0m\u001b[31;49m23.3.1\u001b[0m\u001b[39;49m -> \u001b[0m\u001b[32;49m23.3.2\u001b[0m\n",
      "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m To update, run: \u001b[0m\u001b[32;49mpip install --upgrade pip\u001b[0m\n",
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "%pip install \"pybamm[plot,cite]\" -q    # install PyBaMM if it is not installed\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from numpy import pi\n",
    "\n",
    "import pybamm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "steady-chest",
   "metadata": {},
   "source": [
    "First we will define the parameters in the model. Note the model is posed in non-dimensional form."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "rising-executive",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = pybamm.Parameter(\"Number of winds\")\n",
    "r0 = pybamm.Parameter(\"Inner radius\")\n",
    "eps = (1 - r0) / N  # ratio of sandwich thickness to cell radius\n",
    "delta = pybamm.Parameter(\"Current collector thickness\")\n",
    "delta_p = delta  # assume same thickness\n",
    "delta_n = delta  # assume same thickness\n",
    "l = 1 / 2 - delta_p - delta_n  # active material thickness\n",
    "sigma_p = pybamm.Parameter(\"Positive current collector conductivity\")\n",
    "sigma_n = pybamm.Parameter(\"Negative current collector conductivity\")\n",
    "sigma_a = pybamm.Parameter(\"Active material conductivity\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "worth-easter",
   "metadata": {},
   "source": [
    "Next we define our geometry and model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "virgin-wrestling",
   "metadata": {},
   "outputs": [],
   "source": [
    "# geometry\n",
    "r = pybamm.SpatialVariable(\"radius\", domain=\"cell\", coord_sys=\"cylindrical polar\")\n",
    "geometry = {\"cell\": {r: {\"min\": r0, \"max\": 1}}}\n",
    "\n",
    "# model\n",
    "model = pybamm.BaseModel()\n",
    "phi_p = pybamm.Variable(\"Positive potential\", domain=\"cell\")\n",
    "phi_n = pybamm.Variable(\"Negative potential\", domain=\"cell\")\n",
    "\n",
    "A_p = (2 * sigma_a / eps**4 / l) / (delta_p * sigma_p / 2 / pi**2)\n",
    "A_n = (2 * sigma_a / eps**4 / l) / (delta_n * sigma_n / 2 / pi**2)\n",
    "model.algebraic = {\n",
    "    phi_p: pybamm.div((1 / r**2) * pybamm.grad(phi_p)) + A_p * (phi_n - phi_p),\n",
    "    phi_n: pybamm.div((1 / r**2) * pybamm.grad(phi_n)) - A_n * (phi_n - phi_p),\n",
    "}\n",
    "\n",
    "model.boundary_conditions = {\n",
    "    phi_p: {\n",
    "        \"left\": (0, \"Neumann\"),\n",
    "        \"right\": (1, \"Dirichlet\"),\n",
    "    },\n",
    "    phi_n: {\n",
    "        \"left\": (0, \"Dirichlet\"),\n",
    "        \"right\": (0, \"Neumann\"),\n",
    "    },\n",
    "}\n",
    "\n",
    "model.initial_conditions = {phi_p: 1, phi_n: 0}  # initial guess for solver\n",
    "\n",
    "model.variables = {\"Negative potential\": phi_n, \"Positive potential\": phi_p}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "based-slope",
   "metadata": {},
   "source": [
    "Next we provide values for our parameters, and process our geometry and model, thus replacing the `Parameter` symbols with numerical values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "technological-electric",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pybamm.models.base_model.BaseModel at 0x280b0ed90>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params = pybamm.ParameterValues(\n",
    "    {\n",
    "        \"Number of winds\": 20,\n",
    "        \"Inner radius\": 0.25,\n",
    "        \"Current collector thickness\": 0.05,\n",
    "        \"Positive current collector conductivity\": 5e6,\n",
    "        \"Negative current collector conductivity\": 5e6,\n",
    "        \"Active material conductivity\": 1,\n",
    "    }\n",
    ")\n",
    "params.process_geometry(geometry)\n",
    "params.process_model(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "polyphonic-opinion",
   "metadata": {},
   "source": [
    "We choose to discretise in space using the Finite Volume method on a uniform grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "buried-blind",
   "metadata": {},
   "outputs": [],
   "source": [
    "# mesh\n",
    "submesh_types = {\"cell\": pybamm.Uniform1DSubMesh}\n",
    "var_pts = {r: 100}\n",
    "mesh = pybamm.Mesh(geometry, submesh_types, var_pts)\n",
    "# method\n",
    "spatial_methods = {\"cell\": pybamm.FiniteVolume()}\n",
    "# discretise\n",
    "disc = pybamm.Discretisation(mesh, spatial_methods)\n",
    "disc.process_model(model);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "gothic-deadline",
   "metadata": {},
   "source": [
    "We can now solve the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "straight-anime",
   "metadata": {},
   "outputs": [],
   "source": [
    "# solver\n",
    "solver = pybamm.CasadiAlgebraicSolver()\n",
    "solution = solver.solve(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "excessive-universal",
   "metadata": {},
   "source": [
    "The model gives the homogenised potentials in the negative a positive current collectors. Interestingly, the solid potential has microscale structure, varying linearly in the active material. In order to see this we need to post-process the solution and plot the potential as a function of radial position, being careful to capture the spiral geometry.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "invisible-laser",
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract numerical parameter values\n",
    "# Note: this overrides the definition of the `pybamm.Parameter` objects\n",
    "N = params.evaluate(N)\n",
    "r0 = params.evaluate(r0)\n",
    "eps = params.evaluate(eps)\n",
    "delta = params.evaluate(delta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "affecting-albuquerque",
   "metadata": {},
   "outputs": [],
   "source": [
    "# post-process homogenised potential\n",
    "phi_n = solution[\"Negative potential\"]\n",
    "phi_p = solution[\"Positive potential\"]\n",
    "\n",
    "\n",
    "def alpha(r):\n",
    "    return 2 * (phi_n(r=r) - phi_p(r=r))\n",
    "\n",
    "\n",
    "def phi_am1(r, theta):\n",
    "    # careful here - phi always returns a column vector so we need to add a new axis to r to get the right shape\n",
    "    return alpha(r) * (r[:, np.newaxis] / eps - r0 / eps - delta - theta / 2 / pi) / (\n",
    "        1 - 4 * delta\n",
    "    ) + phi_p(r=r)\n",
    "\n",
    "\n",
    "def phi_am2(r, theta):\n",
    "    # careful here - phi always returns a column vector so we need to add a new axis to r to get the right shape\n",
    "    return alpha(r) * (\n",
    "        r0 / eps + 1 - delta + theta / 2 / pi - r[:, np.newaxis] / eps\n",
    "    ) / (1 - 4 * delta) + phi_p(r=r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "taken-hearing",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define spiral\n",
    "\n",
    "\n",
    "def spiral_pos_inner(t):\n",
    "    return r0 - eps * delta + eps * t / (2 * pi)\n",
    "\n",
    "\n",
    "def spiral_pos_outer(t):\n",
    "    return r0 + eps * delta + eps * t / (2 * pi)\n",
    "\n",
    "\n",
    "def spiral_neg_inner(t):\n",
    "    return r0 - eps * delta + eps / 2 + eps * t / (2 * pi)\n",
    "\n",
    "\n",
    "def spiral_neg_outer(t):\n",
    "    return r0 + eps * delta + eps / 2 + eps * t / (2 * pi)\n",
    "\n",
    "\n",
    "def spiral_am1_inner(t):\n",
    "    return r0 + eps * delta + eps * t / (2 * pi)\n",
    "\n",
    "\n",
    "def spiral_am1_outer(t):\n",
    "    return r0 - eps * delta + eps / 2 + eps * t / (2 * pi)\n",
    "\n",
    "\n",
    "def spiral_am2_inner(t):\n",
    "    return r0 + eps * delta + eps / 2 + eps * t / (2 * pi)\n",
    "\n",
    "\n",
    "def spiral_am2_outer(t):\n",
    "    return r0 - eps * delta + eps + eps * t / (2 * pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "handled-jacksonville",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup fine mesh with nr points per layer\n",
    "nr = 10\n",
    "rr = np.linspace(r0, 1, nr)\n",
    "tt = np.arange(0, (N + 1) * 2 * pi, 2 * pi)\n",
    "# N+1 winds of pos c.c.\n",
    "r_mesh_pos = np.zeros((len(tt), len(rr)))\n",
    "for i in range(len(tt)):\n",
    "    r_mesh_pos[i, :] = np.linspace(spiral_pos_inner(tt[i]), spiral_pos_outer(tt[i]), nr)\n",
    "# N winds of neg, am1, am2\n",
    "r_mesh_neg = np.zeros((len(tt) - 1, len(rr)))\n",
    "r_mesh_am1 = np.zeros((len(tt) - 1, len(rr)))\n",
    "r_mesh_am2 = np.zeros((len(tt) - 1, len(rr)))\n",
    "for i in range(len(tt) - 1):\n",
    "    r_mesh_am2[i, :] = np.linspace(spiral_am2_inner(tt[i]), spiral_am2_outer(tt[i]), nr)\n",
    "    r_mesh_neg[i, :] = np.linspace(spiral_neg_inner(tt[i]), spiral_neg_outer(tt[i]), nr)\n",
    "    r_mesh_am1[i, :] = np.linspace(spiral_am1_inner(tt[i]), spiral_am1_outer(tt[i]), nr)\n",
    "# Combine and sort\n",
    "r_total_mesh = np.vstack((r_mesh_pos, r_mesh_neg, r_mesh_am1, r_mesh_am2))\n",
    "r_total_mesh = np.sort(r_total_mesh, axis=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "monetary-belarus",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArcAAAINCAYAAAAkzFdkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAACCJUlEQVR4nOzdeZxN9R/H8de9s5oZM/Z97LuRsWfLWqJkKQkhqX6KEi2IiAopRZZkTxJCUkSMJFuy78aadcY+Y8w+9/z+OO5EttnvnfF+Ph730dxzz/K5J+nt63O+X4thGAYiIiIiIlmA1dEFiIiIiIikFYVbEREREckyFG5FREREJMtQuBURERGRLEPhVkRERESyDIVbEREREckyFG5FREREJMtQuBURERGRLMPV0QU4ms1m4+zZs2TPnh2LxeLockRERETkPwzD4Nq1axQqVAir9d5jsw98uD179iz+/v6OLkNERERE7uPUqVMUKVLknvs88OE2e/bsgHmzfH19HVyNiIiIiPxXeHg4/v7+ibntXh74cGtvRfD19VW4FREREXFiSWkh1QNlIiIiIpJlKNyKiIiISJahcCsiIiIiWcYD33ObFIZhEB8fT0JCgqNLcRpubm64uLg4ugwRERGRWyjc3kdsbCznzp0jMjLS0aU4FYvFQpEiRfDx8XF0KSIiIiKJFG7vwWazcfz4cVxcXChUqBDu7u5a6AFzJPvChQucPn2aMmXKaARXREREnIbC7T3ExsZis9nw9/fHy8vL0eU4lbx583LixAni4uIUbkVERMRp6IGyJLjfMm8PIo1gi4iIiDNSahMRERGRLEPhVkRERESyDIVbEREREckyFG4fIBUqVGDatGmOLkNEREQk3Wi2hAdEVFQUhw8fpkqVKnfdJzAwkPj4+Nu2//bbbxQqVCg9yxMRERFJEwq3D4i9e/diGAYBAQF33Wfnzp0ZV5CIiIhIOlBbQha3c+dOmjRpQv369bHZbBQtWpSxY8c6uiwRERGRdKGR22QyDHDUSrxeXpCc6WWPHj1Kw4YNeeedd8idOzc2m42aNWvSt29fGjVqRGBgYLrVKiIiIuIITjVyu27dOlq1akWhQoWwWCwsWbLkvsesXbuWatWq4eHhQenSpZk1a1a61hgZCT4+jnklN1T37NmTdu3aMXjwYE6ePEm9evV499138fX15c8//0yfGyQiIiJZ3qJ33uGPceMcXcYdOVW4vX79OlWqVGHixIlJ2v/48eM88cQTNG7cmJ07d/Lmm2/y0ksvsXLlynSu1PmFhISwZs0aevbsSUJCAnv27CEwMBCr1YqLiwvu7u6OLlFEREQyobCTJ+k5ZgyN3nyTX4YOdXQ5t3GqtoQWLVrQokWLJO8/efJkSpQowZgxYwBzqqv169fzxRdf0Lx583Sp0csLIiLS5dRJunZSbd68GZvNRmBgIIcOHSIqKorAwEBOnDjBlStXqFu3bvoVKiIiIlnWqA4duGgYlHN3p/mAAY4u5zZOFW6Ta9OmTTRr1uyWbc2bN+fNN9+86zExMTHExMQkvg8PD0/WNS0W8PZO1iEOERsbC0B0dDQ7duygWLFi5MqVi9GjRxMQEEDlypUdXKGIiIhkNqf++ouxmzcDMPrtt3HLls3BFd0uU4fbkJAQ8ufPf8u2/PnzEx4eTlRUFNnucMNHjhzJsGHDMqpEh6lTpw6urq4MHz6ciIgISpYsyYQJExg/fjzr1q1zdHkiIiKSCQ3u1Ilo4BE/P1p9+KGjy7mjTB1uU2LgwIH069cv8X14eDj+/v4OrCh9+Pv7M2PGDPr378+5c+dwdXUlMjKSFStWUL16dUeXJyIiIpnMtx9+yOxjxwD4bNw4LFanenQrUaYOtwUKFCA0NPSWbaGhofj6+t5x1BbAw8MDDw+PjCjP4bp06UKXLl3IlSsXs2bN4qmnnnJ0SSIiIpIJGTYbbw4fDkDt3Lmp2a2bgyu6O+eM3ElUp04dgoKCbtm2atUq6tSp46CKnM/p06e5cuXKPVcmExEREbmX+GXLMOLjAegzZIiDq7k3pwq3ERER7Ny5M3EZ2OPHj7Nz505OnjwJmC0FXbt2Tdy/Z8+eHDt2jHfffZeDBw8yadIkFixYQN++fR1RvlPas2cP3t7elChRwtGliIiISGYUH4/LgAFcv/G2jpP/TbBThdutW7dStWpVqlatCkC/fv2oWrUqQ278CeHcuXOJQRegRIkSLFu2jFWrVlGlShXGjBnDtGnT0m0asMyoRYsWREREYEnO0mYiIiIidrNmEbJ/P7GAi4sLRYoUcXRF9+RUPbeNGjXCMIy7fn6n1ccaNWrEjh070rEqERERkQfU9eswZAgnbrwtUqQIrq5OFR9v41QjtyIiIiLiRMaMgXPnOJ43L0CmaHNUuBURERGR24WEwOjRAJxo0gSA4sWLO7CgpFG4FREREZHbDR1qtiXUrs3xG8uzauRWRERERDKf/fth2jTz5zFjOPHPP4BGbkVEREQkM3r3XbDZoG1bqFeP48ePAxq5FREREZHM5vffYdkycHWFUaNISEhInIq1WLFiDi7u/hRuRURERMSUkABvv23+3LMnlC3L6dOniY+Px93dncKFCzu2viRQuBURERER08yZsH07+PrCjUW0jh07Bpj9ti4uLo6sLkkUbkVEREQErl6FgQPNn4cNgxtz29rDbWbotwWFW7mDWbNm3XE1OBEREcnChg6FixehYkXo1Stxsz3clixZ0lGVJYvC7QOkQoUKTLNP6yEiIiJit3cvTJxo/jxuHLi5JX6kcCtOKSoqisOHD1OlShVHlyIiIiLOxDCgTx/zYbJ27aBZs1s+zmzh1tXRBUjG2Lt3L4ZhEBAQcMfPY2NjqVWrFgCXL18GYOzYsQBs2bIFd3f3DKlTREREMtiiRbBmDXh6wpgxt32scCtOZefOnfTr148NGzZgs9koWrQogwYN4s0337xlP3d3d3bu3AmQ2G/7wgsvZGitIiIiksEiI+Gtt8yf+/eH/6xAFh4ezsWLF4HM80CZwm0WdvToURo2bMg777xD7ty5sdls1KxZk759+9KoUSMCAwMdXaKIiIg40iefwMmTULSouSrZf9hXJsudOzd+fn4ZXV2KKNwml2GYf8pxBC8vsFiSvHvPnj1p164dgwcPpnbt2nTo0IE333yTkSNH8ueffyrcioiIPMiOHDHDLZjtCF5et+1y9OhRIPO0JIDCbfJFRoKPj2OuHREB3t5J2jUkJIQ1a9awceNGEhIS2LNnDyNHjsRqteLi4nLPHlq1I4iIiGRxhgGvvgoxMfDoo/D003fc7ciRIwCULl06I6tLFc2WkEVt3rwZm81GYGAghw4dIioqisDAQE6cOMGVK1eoW7euo0sUERERR5k7F1avNh8i++qru/7NsH3kNjOFW43cJpeXlzmC6qhrJ1FsbCwA0dHR7Nixg2LFipErVy5Gjx5NQEAAlStXTq8qRURExJldvgx9+5o/v/8+lCp1110z48itwm1yWSxJbg1wpDp16uDq6srw4cOJiIigZMmSTJgwgfHjx7Nu3TpHlyciIiKOMmAAXLhgrkT29tv33NU+clvqHgHY2SjcZlH+/v7MmDGD/v37c+7cOVxdXYmMjGTFihVUr17d0eWJiIiII6xfD1Onmj9//TXc4xmcmJgYTp48CWSukVv13GZhXbp04ezZs+TMmZNFixaxefNmGjRo4OiyRERExBFiY+F//zN/fuklqF//nrsfP34cwzDw8fEhX758GVBg2lC4zeJOnz7NlStX7roymYiIiDwgPvsM9u+HvHn/nQLsHm5+mMySjKlIHU3hNovbs2cP3t7emWZVEREREUkHhw/Dhx+aP3/xBeTKdd9D7A+TZaZ+W1DPbZbXokULIhw1u4OIiIg4ns0GL74I0dHQrBl06pSkww4fPgxkrn5b0MitiIiISNY2YYL5IJmPj/kwWRJbDIKDgwEoV65celaX5hRuRURERLKqI0fMqb8APv0UihdP8qH2cFumTJl0KCz9KNyKiIiIZEU2G/ToAVFR0KQJvPJKkg+Njo5OnAasbNmy6VVhulC4FREREcmKJk2CdevMxaemTQNr0mPf0aNHMQwDPz8/8ubNm45Fpj2FWxEREZGs5tgx6N/f/Hn0aEjmrEn2loSyZctmqmnAQOFWREREJGuxtyNERkKjRtCzZ7JPkVn7bUHhVkRERCRr+eorWLvWbEeYPj1Z7Qh29mnAMlu/LSjcioiIiGQdx4//247wySdQsmSKTnPo0CFAI7ciIiIi4igJCdC1K1y/Dg0bwquvpvhU9nBboUKFtKouwyjcioiIiGQFo0aZizVkzw4zZ6aoHQHg0qVLXLhwAVBbgoiIiIg4wP7581k/ZIj5ZsKEZM+OcDP7qK2/vz/e3t5pUV6GUrgVERERycT2bdlC4HPP0dBmY33jxtClS6rOd/DgQQDKly+fFuVlOIXbB0iFChWYNm2ao8sQERGRNDThhReIA2zAiqpVIZXz0ircSqYQFRXF4cOHqVKliqNLERERkTTy86BBTD5wIPH9pJkzCQsLS9U5FW4lU9i7dy+GYRAQEODoUkRERCQNrFu0iBdGjACgT7VqlC9fnitXrjBp0qRUnVfhVpzazp07adKkCfXr18dms1G0aFHGjh3r6LJEREQkFWyxsbR97jkuA/5ubnzy++8MHDgQgLFjxxIVFZWi88bExHDs2DEAypUrl1blZihXRxeQ2RiGQWRkpEOu7eXllaz1nY8ePUrDhg155513yJ07NzabjZo1a9K3b18aNWpEYGBg+hUrIiIi6eb6kCGEx8cDMPbzz/Hw9aVjx44MGTKEf/75h1mzZvFqCua5DQ4OJiEhAT8/PwoVKpTWZWcIjdwmU2RkJD4+Pg55JTdU9+zZk3bt2jF48GBOnjxJvXr1ePfdd/H19eXPP/9MpzskIiIi6WrNGo5+8gnxgLeHB61vhFg3NzfeeustAL788ksMw0j2qffv3w9AxYoVkzWg5kwUbrOokJAQ1qxZQ8+ePUlISGDPnj0EBgZitVpxcXHB3d3d0SWKiIhIcp0/D507s+HG23oNG+Li4pL48QsvvED27Nk5ePAgq1evTvbp9+3bB5jhNrNSW0IyeXl5ERER4bBrJ9XmzZux2WwEBgZy6NAhoqKiCAwM5MSJE1y5coW6deumY6UiIiKS5mw2c3ndkBA2+PpCeDj16tW7ZZfs2bPzwgsvMH78eMaPH8+jjz6arEvYR24rVaqUZmVnNIXbZLJYLJlitY7Y2FgAoqOj2bFjB8WKFSNXrlyMHj2agIAAKleu7OAKRUREJFk++wxWrgRPTzZ4e98x3AL06tWL8ePHs2zZMs6cOUPhwoWTfImb2xIyK7UlZFF16tTB1dWV4cOHs27dOkqWLMmECRMYP348s2bNcnR5IiIikhybN8OgQQCc+uADTp47h4uLC7Vr175t13LlyiXOkvTNN98k+RKxsbEcPnwYyNzhViO3WZS/vz8zZsygf//+nDt3DldXVyIjI1mxYgXVq1d3dHkiIiKSVBcvwrPPQnw8PPccf9yYxaB69er4+Pjc8ZAePXqwfv16ZsyYwcCBA5P0cNjhw4eJj48ne/bsFClSJE2/QkbSyG0W1qVLF86ePUvOnDlZtGgRmzdvpkGDBo4uS0RERJIqIQE6dYJTp6BMGfj6a/5Ytw6ARx555K6HtW/fHh8fH44ePcrGjRuTdKm9e/cCmXumBFC4zfJOnz7NlStXtDKZiIhIZjRsGKxaBV5esHgx+Pqy7ka4bdiw4V0P8/b2pl27dgB89913SbqUPdxm9udyFG6zuD179uDt7U2JEiUcXYqIiIgkx7Jl8OGH5s9TpkBAAOfOnSM4OBiLxXLHh8lu1rlzZwAWLFhAXFzcfS+3Z88eQOFWnFyLFi2IiIjI1H+9ICIi8sA5fhyef978uVcvuBFU//jjDwAeeughcubMec9TNGnShPz583Pp0iV+++23+17SPnKb2f+2V+FWRERExJlERcHTT8PVq1C7Nnz+eeJHv//+OwCNGze+72lcXV1p3749AAsXLrznvtevX+fYsWOARm5FREREJK0YhjlSu2MH5MkDP/wAN60qag+3TZo0SdLpnnnmGQB++umne7Ym7N+/H8MwyJcvH3nz5k3FF3A8hVsRERERZzFxIsycCVYrfP89+PsnfnT69GkOHz6M1Wq950wJN6tfvz758uXjypUricH4TrJKvy0o3CaJYRiOLsHp6J6IiIiksbVr4c03zZ9Hj4ZmzW752B5Oq1evjp+fX5JO6eLiQtu2bQFYvHjxXffbvXs3oHCb5bm5uQEQGRnp4Eqcj315XxcXFwdXIiIikgWcOAHt25vz2j7/PPTrd9suQUFBQNL6bW/WunVrAH755Ze7Dk7t2rULgCpVqiTr3M5IK5Tdg4uLCzly5OD8+fMAeHl5adYBwGazceHCBby8vHB11S8hERGRVLl+Hdq0MVciq17dnPbrP3nDMAxWr14NQLP/jOjeT+PGjfH29ubMmTPs2LGDatWq3XZuhdsHSIECBQASA66YrFYrRYsWVdgXERFJDcOAHj1g1y7Ilw9+/BGyZbttt0OHDnHmzBk8PDyoX79+si7h6enJo48+ypIlS/jll19uC7enTp3iypUruLq6UrFixVR9HWegcHsfFouFggULki9fviRNgPygcHd3x2pVV4uIiEiqfPIJzJ8Prq6wcOEtD5DdbNWqVYD5gFi2O4Tf+3nyySdZsmQJP//8M0OGDLnlM/uobYUKFfDw8Ej2uZ2Nwm0Subi4qL9URERE0s7SpfDee+bP48dDgwZ33dXekvDoo4+m6FItW7YEYNu2bVy8eJE8efIkfpaVWhJAD5SJiIiIZLxdu6BTJ7MtoWdP83UXsbGxrFmzBkh5uC1YsCCVK1fGMIzEUWC7nTt3Agq3IiIiIpISISHQqpX5IFnTpvDll/fcfePGjURERJAvXz4CAwNTfNnmzZsDsHLlylu279ixAyBV53YmCrciIiIiGSUqypwZ4dQpKFvWXIHsxtSjd7NixQrADKeped7FHm5/++23xCnBrly5krjs7n8fNMusFG5FREREMoJhwIsvwl9/Qc6c8Msv5j/vwx5uH3/88VRd3v4w2rlz5zhw4ADwb0tC8eLFyZUrV5LOYxgQGgobN5o/Oxs9UCYiIiKSEYYPh3nzzJkRFi+GMmXue8iZM2fYtWsXFoslxf22dp6entSvX59Vq1axZs0aKlasyPbt24HbR21jY+Gff+DoUfN17Ni//zx2zOyoADh3Dm7Mmuo0FG5FRERE0tucOfDBB+bPX30FjRol6bDly5cDUKtWLfLmzZvqMpo0aZIYbnv37p0YbosVq87MmbBhg/kKDgab7e7nsVjMWcsuXlS4FREREXmwrFljtiMAvP02vPRSkg9dtmwZAE888USalGJfujcoaA1duoxjwYKlAHzxxe0zJXh5QcmSUKrU7f8sVgycdUpchVsRERGR9LJ3L7RtC3Fx8Oyz5qINSRQdHZ04bdeTTz6Z6lIuX4bdu6tisbgQHh7GnDlvJn7m6mrl4YehXj3zVaOGOSKbGRciVbgVERERSQ9nzkCLFhAeDvXrwzffQDJmO1i7di2RkZEUKlQoxdN0xcWZK/rOmAFBQRAf7w7kAUJxdS1FYGALAgJyMn58M3x8UnQJp6NwKyIiIpLGws+cwbtlS1xOn4Zy5eCnn8DTM1nnWLrUbBl48sknsSRzCPXMGZgyBaZONR/6snvoIahceQatWuWlQ4eayTpnZqFwKyIiIpKG4qOiKFWiBNni4liRKxcVly+HJE6zZWcYRmK4bd26dZKP27wZxowxR2sTEsxt+fPDyy9Dly7m1LrQMlm1ZDZON8/txIkTKV68OJ6entSuXZstW7bcc/+xY8dSrlw5smXLhr+/P3379iU6OjqDqhURERH5l2Gz0bFiRS7GxXEKqBIWxnvTphEVFZWs82zbto0zZ87g7e1NkyZN7n1NA1avNhc7q1MHFi40g22DBubMYydPwocf2oNt1udU4Xb+/Pn069ePoUOHsn37dqpUqULz5s05f/78HfefO3cuAwYMYOjQoRw4cIDp06czf/583nvvvQyuXERERARGNm/OwhMnAKhSvDjxCQmMHDmS6tWrs2vXriSf56effgLMhRs879LOYBhmt0Pt2vDoo+akDK6u8MILsHs3rFsHHTqAu3tqv1Xm4lTh9vPPP+fll1+me/fuVKxYkcmTJ+Pl5cWMGTPuuP/GjRupV68enTp1onjx4jz22GN07NjxvqO9IiIiImmtbbVqDFq9GoAv27dnx7FjLFmyhIIFC3LgwAHq1q3LokWLknQu+35t2rS54+erV5uhtk0b+PtvyJYNXn/dXGhh5kyoXDktvlHm5DThNjY2lm3bttGsWbPEbVarlWbNmrFp06Y7HlO3bl22bduWGGaPHTvG8uXLadny7r0kMTExhIeH3/ISERERSY1zM2bw044dALQrW5bXFyzAYrHQunVrdu/ezWOPPUZkZCTt27fnq6++uue5Dhw4wIEDB3Bzc7ttCrC//jLbDx591Ay13t4wcCCcOAFffglFi6bXN8w8nCbcXrx4kYSEBPLnz3/L9vz58xMSEnLHYzp16sTw4cOpX78+bm5ulCpVikaNGt2zLWHkyJH4+fklvvz9/dP0e4iIiMgD5o8/WP+//2EA3q6ufL979y0f58mTh2XLlvHaa69hGAavvfYaEyZMuOvpFi9eDEDTpk3JkSMHAAcOmKO0Dz9sth+4u8Mbb5gjtSNGQL586fTdMiGnCbcpsXbtWkaMGMGkSZPYvn07ixcvZtmyZXz44Yd3PWbgwIGEhYUlvk6dOpWBFYuIiEiWsmMHPPUU0+PjAejz9tu432HpLldXVyZMmMCAAQMAeP3115kzZ84dT7lw4UIAnn76aS5fNkNs5cpmf63VCt27m8vjjhtnzoQgt3KaqcDy5MmDi4sLoaGht2wPDQ2lwF0WLX7//ffp0qULL91Yxq5y5cpcv36dV155hUGDBmG9w0TJHh4eeDjrenEiIiKSeRw5Ao8/zj/h4fx2Y9OL91ha12KxMGLECKKjoxk7diwvvvgiRYsW5ZFHHkncJzg4mJ07d+Li4sKlS20oU8ZcWQygdWsYORIqVEjH75QFOM3Irbu7O9WrVycoKChxm81mIygoiDp16tzxmMjIyNsCrIuLC2DODyciIiKSLs6dg+bN4fx5pubLh4HZRlCqVKl7HmaxWBgzZgxPP/00cXFxPPPMM5w5cybx8wULFgCQLVszBgzIw+XLEBBgPkC2ZImCbVI4TbgF6NevH1OnTuWbb77hwIEDvPrqq1y/fp3u3bsD0LVrVwYOHJi4f6tWrfjqq6+YN28ex48fZ9WqVbz//vu0atUqMeSKiIiIpKmrV+Hxx+HYMWJLlmTajQG1V199NUmHW61WZs+eTZUqVbhw4QIdOnQgPj6erVtPMnz4pwBERDxD7twwaZLZ+dC0aXp9mazHadoSADp06MCFCxcYMmQIISEhBAYGsmLFisSHzE6ePHnLSO3gwYOxWCwMHjyYM2fOkDdvXlq1asXHH3/sqK8gIiIiWVlUFDz1lDmRbIECLHnzTULfeIMCBQrw1FNPJfk0Xl5eLFy4kOrVq7NhwwYqV36CgwfXAeZCVE8+WZjZsyFnznT6HlmYxXjA//4+PDwcPz8/wsLC8PX1dXQ5IiIi4qzi4+Hpp2HpUvDzgz/+oOEbb7Bu3ToGDx58zwfa72bw4Kl8/PErie+9vGrw4otPM378gLSsPNNLTl5zqpFbEREREaeUkADdupnB1tMTfv6Z3RYL69atw9XVlZ49eybrdFFRMHQofPbZS8CHwGkaN+7Db799jqurJV2+woNC4VZERETkXmw2eOUVmDvXXN/2hx+gQQPGv/wyAO3ataNw4cJJPt369fDii3D4MICFJ55YTd++UTRtWiV96n/AKNyKiIiI3I1hmBPNzphhTjL7/ffw5JNcvHgxcZ7a119/PUmnioiA996DCRPM0xYqBJMnQ6tWZdPzGzxwFG5FRERE7sQw4N13YeJEsFjgm2/gmWcAmDp1KtHR0VSvXp169erd91SrV8PLL5vL5AL06AGffQY3FiCTNKRwKyIiInInZlOs+fOUKfD88wDExsYmLp/bp08fLJa798iGhcHbb8O0aeb7YsVg6lR49NF0rfyB5lTz3IqIiIg4hZEjwT77wZdfwk0rj33//fecPXuWQoUK0aFDh7ue4tdfoVKlf4Nt796wd6+CbXrTyK2IiIjIzcaONZtjAT75BG7qqTUMg89ujOb26dMHd3f32w6PjIR33jEXYAAoUwamT4cGDdK7cAGFWxEREZF/ff019O1r/vzBB2bP7U1+/fVX9u7di4+PD6+88spth2/bBp07w6FD5vs334QRIyBbtvQtW/6ltgQRERERMB8Ys89X278/DBly2y6jRo0CoGfPnuS46WmwhASzk+Hhh81gW6gQrFoFX3yhYJvRNHIrIiIiMn++OfksmFN/jRxpzpBwkw0bNvDnn3/i7u5OX/voLuYMCF26mPPXgjmhwtdfQ65cGVS73EIjtyIiIvJgW7LEnAnBZjPn6xo79rZgCyQur9utWzcKFSqEYcDs2fDQQ2awzZ7dHPxdsEDB1pE0cisiIiIPrkWL4LnnID7eDLiTJ98x2G7ZsoWVK1fi4uLCgAEDuHzZ7GD44Qfz83r14NtvoUSJDK5fbqORWxEREXkw/fADdOhgBtvOnWHmTHMVsjsYPnw4AJ07d+bYsZJUrmwe7uoKH38Mf/yhYOssNHIrIiIiD55588yR2oQE6NrVXF7XxeWOu27ZsoVly5bh4uKCh8fgxHlqy5aF776DGjUysG65L43cioiIyIPlu+/MkdqEBHjhhXsGW4ChQ4cCkDPn80ydWgaA//0Ptm9XsHVGGrkVERGRB8fs2WagNQxz1bGvv75rKwLAn3/+yYoVKwBXLl58nxw5zBXHnn46owqW5FK4FRERkQfDzJnQo4cZbP/3P3MJsXsE2/PnI2jevP2Ndy9St24p5s6FYsUyplxJGYVbERERydIiIiIY9vTTXP7tN3yAkFKluPbPP9hatsTT05NcuXLh7+9P+fLlqVGjBqVLl2bx4j0899xTxMeHAtC+fSvmzjUfIBPnpn9FIiIikqXN7NGDz3777d8NR4+ar7vIli0HUVHXgAQslty0a9eVBQueTP9CJU0o3IqIiEiWZBgGs7p25Y0FCwAo4evLM//7H4ULF8bX1xcXFxeioqK4ePEi//zzDzt37mHbtu1ERV0FwM+vNps2/UKFCnkc+C0kuRRuRUREJEt6MiCA5fv3A/BqlSpM2LYN611mRdi7F9q3B5stGotlNMWLr2b//tV4erpnZMmSBjQVmIiIiGQ5f735ZmKwfaJkSSZu337XYDtrFtSqBQcPQuHCnqxbN4Rjx9Yp2GZSCrciIiKStYwcySfjxgFQ3M+PpcHBWO4wK0JkJHTvbr6ioqB5c9ixA+rXz+iCJS0p3IqIiEjWYBgwYAC/vfcePwIuFgu/rF9/xxHbAwfM0dpZs8zZwD76CJYvh7x5M7xqSWPquRUREZHMz2aD3r2J+eoret/Y1PuNN6gUEHDbrnPmQM+ecP06FCgA338PjRplaLWSjhRuRUREJHOLj4cXX4Rvv+VT4DBQoEABhg0bdstuUVHwxhvmCmMATZuaK/Hmz5/hFUs6UluCiIiIZF4xMfDss/DttxyxWvnIzQ2AMWPG4Ofnl7hbcDA8/LAZbC0W+OADWLlSwTYrUrgVERGRzOn6dXjqKfjxRww3N16rXJmYuDiaNWtGx44dE3dbuBCqV4fduyFfPli1CoYOhbtMniCZnMKtiIiIZD5hYeb0Br/9Bl5ezOnbl1W7duHh4cGkSZOwWCzExcFbb5nz10ZEQMOGsHOn2Y4gWZfCrYiIiGQuFy5AkyawYQP4+XF+wQL6Tp8OwNChQylTpgznzpkh9vPPzUP694fVq6FgQQfWLRlCD5SJiIhI5nH2LDRrZs7llTcv/PYbb4waxaVLl6hSpQpvv/02f/5ptuGGhICvL3zzDbRp4+jCJaNo5FZEREQyh6NHzRUWDhyAIkVg3ToWHzvG/PnzcXFxYfr0GUyY4EbjxmawDQiAv/9WsH3QaORWREREnN+OHfD443D+PJQsCUFBXPTx4dVXXwXgzTffZfToaixYYO7eqRNMmQLe3g6sWRxC4VZERESc2++/Q+vWcO0aBAbCr79i5M9Pz/btOX/+PKVLV+KXX4Zy6BC4usIXX0CvXuaUX/LgUbgVERER57VokTkMGxtrLiO2ZAn4+fHdnDksWrQIFxdXzpyZTVSUB4ULww8/QJ06ji5aHEk9tyIiIuKcvv7anMcrNhbatYNffwU/P/755x969eoFQELCEKKiqtG4MWzfrmArCrciIiLibAwDhg+Hnj3Nn195BRYsAE9PEhIS6NChC+Hh4cDDwED69zenu82Xz9GFizNQW4KIiIg4j4QE6NMHJk4037//PgwblthA26JFN/7660/ABx+fb/n2W1fNhiC3ULgVERER5xATA127mqO0Fgt8+SX07g1AfHwCdeu+xN9/fwdA/vyDWbeuNGXLOrJgcUYKtyIiIuJ4166ZfbWrV4ObG3z7LXToAMCuXUdo1OhZrl7dAUD27GXZv/9NcuVyZMHirNRzKyIiIo5lX0539Wrw9iZ+6VL+KFCAd999l9KlKxEYWOZGsPWka9c5XLy4l1y5PBxdtTgpjdyKiIiI45w4Ac2bYwQH86efH3ObNGFRly5cvHjxPzu68f77Uxg+vLMjqpRMROFWREREHGP7di60aMGM8+eZ6urK0bAw+PFHADw9cxEd/QTQgtq1azBnTj5Kl/ZzbL2SKSjcioiISIab+sYbLJk0iaCEBGIA4uPx8fGhZctn2LevI/v2NQbc6N8fPvrIXHlMJCn0S0VEREQy1MCWLRn166+J72tUrcprr79O4cLP0q2bNyEh4OsL33yDpvmSZNMDZSIiIpIhoqOi+Khx48Rg6+/lxdpVq9iybTtXrnSnZUsz2AYEwN9/K9hKymjkVkRERNJdQnQ0VQsW5GBYGAB9atXi840buR7pQocO8MMP5n6dOsGUKeDt7cBiJVPTyK2IiIikr2vX+PyhhxKD7WsNGzL2r784FOxCrVpmsHVzgwkTYM4cBVtJHY3cioiISPo5d44tDRsy+PBhADo1bMjEtWtZsABefBGuX4fChc2AW6eOg2uVLEEjtyIiIpI+9u8ntGZN2h0+TCzQumFDZqxcQ9++5uJj169D48awfbuCraQdjdyKiIhI2lu3jpinnqJdWBhngPKlSjF6ylKaNbOyfr25i6b5kvSgX04iIiKStubMwejenV7x8WwE/Hx9GTxsOY884ktoqKb5kvSlcCsiIiJpwzBg2DAYNoyxwHTAarXS/tn5dOtWloQEc5qvxYuhTBkH1ypZlsKtiIiIpF5MDPToAd99x3LgbYsFDIOAgM+YNu1xADp3hq+/1mwIkr4UbkVERCR1Ll2Ctm3hzz/Z7eLCc25u2KKj8fPrwe7db+LmBl98Aa+9BhaLo4uVrE7hVkRERFIuOBieeAKOHOFs9uw86enJtQsXsFobERY2icKFLZrmSzKUpgITERGRlFm3zkytR45wrWhRnihchFMXLgDlsNkW07ixu6b5kgyncCsiIiLJ9+230KwZXL5MXM2atPYvxc6DB4C8wHL698/Jb79BvnyOLlQeNGpLEBERkaQzDPjgAxg+3Hz7zDO0uZaN31d+C3jh7f0Lc+aU1DRf4jAKtyIiIpI00dHmjAhz5wJg6z+ACou3E3x4IWClaNH5rFpVi7JlHVumPNgUbkVEROT+QkLMGRE2bwZXVy59Mo7yIxdw8eIfAFSp8hEbNjypab7E4RRuRURE5N527oSnnoJTpyBnTua9NIgu/UcQH38GsFChQnO2bx+AVU/yiBPQL0MRERG5uyVLoF49OHWKI8WLU79oLTp++jbx8WewWHIzYcLv7Nu3HKtVE9iKc9DIrYiIiNzOMGDUKHjvPU4AwwoW5Jt/TmEYJwBwcfFj/vwfePrpho6sUuQ2CrciIiJyq+hoePllQufM4SPga6uVuHPnbnzYgueee5cZM+qRLZubI6sUuSOFWxEREflXaCgRrVox5u+/+RS4DmCzAc3w8vqQ7757WNN8iVNTuBUREREAErZv55tmzRh05QohiVtrASMJCGjCokVomi9xenqgTERERBjZpQsVatSgx41g6+dVFFgAbKZTpyZs3qxgK5mDRm5FREQeYIbNRtfq1ZmzcycA2V1cyJFzKKcu9sfV1Z0vvoBevcCiyRAkk1C4FREReUCdOXyY1xo1YunZswCU8c5JiG0Hpy4Wo1AhWLgQ6tRxcJEiyaS2BBERkQdQ8O+/U7ZcOZaePYsL8EKZZzl8/TLXoorRqBFs365gK5mTwq2IiMgD5vqKFXR/7DEiDQML0KDQAGYdng/Au+/CqlWQP79jaxRJKacLtxMnTqR48eJ4enpSu3ZttmzZcs/9r169Sq9evShYsCAeHh6ULVuW5cuXZ1C1IiIimYhhEDVuHK1btmRjfDw+FivFfKax9uxIfH1h8WL45BNwVdOiZGJO9ct3/vz59OvXj8mTJ1O7dm3Gjh1L8+bNOXToEPny5btt/9jYWB599FHy5cvHwoULKVy4MP/88w85cuTI+OJFREScWWws0a++SrsZMwgCvKyuxNl+40REY6pUMftrS5d2dJEiqWcxDMNwdBF2tWvXpmbNmkyYMAEAm82Gv78/r7/+OgMGDLht/8mTJ/Ppp59y8OBB3NxStkpKeHg4fn5+hIWF4evrm6r6RUREnFJoKDFt29Ju0yaWA+4WN2KNVUBDevSA8eMhWzZHFylyd8nJa07TlhAbG8u2bdto1qxZ4jar1UqzZs3YtGnTHY9ZunQpderUoVevXuTPn5+AgABGjBhBQkJCRpUtIiLi3LZtI6Z6dZ6+EWxd8CDWWImnZ0NmzIBp0xRsJWtxmraEixcvkpCQQP7/dLDnz5+fgwcP3vGYY8eOsWbNGjp37szy5cs5cuQIr732GnFxcQwdOvSOx8TExBATE5P4Pjw8PO2+hIiIiDOZO5eYF1/k6ZgYlgEWPEhgOaVLN2bhQqhSxdEFiqQ9pxm5TQmbzUa+fPmYMmUK1atXp0OHDgwaNIjJkyff9ZiRI0fi5+eX+PL398/AikVERDJAfDy8+y7RnTvT9kawBU8MltGuXRO2blWwlazLacJtnjx5cHFxITQ09JbtoaGhFChQ4I7HFCxYkLJly+Li4pK4rUKFCoSEhBAbG3vHYwYOHEhYWFji69SpU2n3JURERBzt4kV4/HEiP/2Up4BfAciG1bqMzz9vysKF4Ofn2BJF0pPThFt3d3eqV69OUFBQ4jabzUZQUBB17jKLdL169Thy5Ag2my1xW3BwMAULFsTd3f2Ox3h4eODr63vLS0REJEvYtg2qVyciKIiWFiurAPAmd+7lrFvXhL59tYyuZH1OE24B+vXrx9SpU/nmm284cOAAr776KtevX6d79+4AdO3alYEDBybu/+qrr3L58mX69OlDcHAwy5YtY8SIEfTq1ctRX0FERMQxZs2CevW4cPIkjVw8+MOwAdmpXn0l+/c3ol49RxcokjGc5oEygA4dOnDhwgWGDBlCSEgIgYGBrFixIvEhs5MnT2K1/pvH/f39WblyJX379uWhhx6icOHC9OnTh/79+zvqK4iIiGSs2Fjo2xcmTWI3UBsXohNigBy8+OJKpkypxU3deyJZnlPNc+sImudWREQyrbNnoX172LiRFUAb3IkhFrDw4Ye/Mnhwc0dXKJImkpPXnGrkVkRERJJo/Xpo356wkBBaWFzYZCQAsbi7F2f69Ck8//yjjq5QxCEUbkVERDITw4BJkzD69GFBQgKvW1y5YMQDUKxYS3bvXoSvr6eDixRxHKd6oExERETuISoKundnb+/eNE5I4Dm4EWxz06LFmxw//ouCrTzwNHIrIiKSGRw/TnjbtnywaxdfAuZC857kyzeQlSvfITBQa+iKgEZuRUREnJ7xyy/MrVyZcrt28QX2YNuOZ589yPHjQxRsRW6icCsiIuKsEhI4+L//0bRVKzpfv04I4EoJPDxWMHv2IubPL4aXl6OLFHEuaksQERFxQmd27WJsq1aMO3WKOMANFxIYQpkK/Vm0yIMKFRxdoYhz0sitiIiIk5kxcCDFAwP57EawLURV4jjMCy8OYetWBVuRe9HIrYiIiJNIiI/nqw4d6Ld4MfGAC5CPD7mSbRDfTLbQtaujKxRxfgq3IiIiTiDs1CmqVajAsevXAShNLs6ygZwVy7P6B6hY0cEFimQSaksQERFxsOtbtvBUuXKJwbYidTjCeZ59oTxbtijYiiSHwq2IiIgDXfrqK5rWqcO6qCjcAX9e43i2jcyc6cLMmeDt7egKRTIXtSWIiIg4QnQ0//ToQfO5czkEZMcNK0vwqdCSX3+ASpUcXaBI5qRwKyIiktFOnGB3y5Y8fuAA5wBffLnGBrp0DWDSJI3WiqSG2hJEREQy0s8/83vlyjS4EWx9KEas5z5mzAzgm28UbEVSS+FWREQkI8TGwltvMe+pp2geEUE44EEtCpbZyd9/F+GFFxxdoEjWoHArIiKS3k6cwGjQgE8//5yOQBxgoS1Pd/qD7dtzEBDg6AJFsg6FWxERkfS0ZAkJgYG8sWUL797Y5Orahxkzf+C77zzx8XFodSJZjh4oExERSQ+xsfDuu4SPG8fzwM8AWChUaAy//96XsmUdW55IVqVwKyIiktaOH4cOHdj+9980BCIA8KBFizn8+OMzeHg4tjyRrExtCSIiImlp8WKoWpXFf/9NA+zB1sKgQQtZvlzBViS9aeRWREQkLcTEwNtvEzphAs2BXTc2u7sXZebMGXTq1NSR1Yk8MBRuRUREUuvoUYxnn2Xu9u30AS7d2FypUge2bJmNl5e7I6sTeaCoLUFERCQ1fviBfVWq0Hj7dp7HDLZWa35efPEj9u6dp2ArksE0cisiIpIS168T0bs3w2fN4gsgHrDgSZmy7xMU9BZFiqi5VsQRNHIrIiKSTMbOnSwsW5YKs2bxKWawtVpaMWjwAQ4ceE/BVsSBNHIrIiKSVIbBkSFD6P3xx6w0DAA8KUiOAlP4+ecnqVHDwfWJiMKtiIhIUpzbu5dxbdrwxdGjxAKuWHGjH0+2H860adnw9XV0hSICCrciIiL3NfGNN+gzfjwJN977U54LHj8x8auyvPACWCyOrE5EbqZwKyIichdR4eF82Lw5n2zejA1wAQownByVB7NyvoUKFRxdoYj8l8KtiIjIHexZtowGrVsTlmCO15amAKdYT/s3SzFyJHh6OrhAEbmjVIXb3bt38+eff+Lu7k7dunWpVKlSWtUlIiLiMAc//ZRW/fsTduOhsQo8y6V881nyDTz+uIOLE5F7SnG4HTduHH379sXX1xcXFxeuXLlC5cqV+eabbwgMDEzDEkVERDLItWusbteOZ1avJgzIhRsWplC8xQv8PhPy53d0gSJyP8ma53bGjBls376dmJgYPv74Y0aNGsWVK1e4dOkSx44do0WLFjRo0ICNGzemV70iIiLpY+tWJpUqxeM3gq0//kS5nWDoly+wbJmCrUhmYTGMG3/nkgQVK1bk8OHDANhsNtq1a0e9evWoWrUqgYGB+Pn5MXHiRL777rtME3DDw8Px8/MjLCwMX83jIiLy4LHZiP/0U/q89x6TbDYA8vI4eSosYf58DypXdnB9IpKsvJaskdv9+/dz7do1Nm7ciJubG1arlXnz5tGyZUty5cpFyZIl+fHHH9m2bRvLli3jxIkTqfkeIiIi6ev0aa42akTLAQOYZLNhATwZwrO9lrNtm4KtSGaU7J5bT09PatasSb169ahSpQrz58/HZrNx8OBBdu7cybp161izZg1du3blypUr+Pj4EB4enh61i4iIpNyCBQS/9BJPXbvGIcAVd7L5fM/c79vx5JOOLk5EUirFD5SNGTOGRo0acezYMXr27EmVKlXw9/dn+/btFCpUiNOnT3P69Gn27t2blvWKiIikTlgY9O5N0Jw5tAeuAK4UoMbDv7J4cSAFCzq6QBFJjRSH28DAQLZt20bPnj15+OGHsbfuurq6MmPGDACKFClCkSJF0qZSERGR1Fq3Drp04auTJ3kdSAAs1GLQ0J8YMqQA1mQ164mIM0rVPLelSpVi1apVhIaGsnnzZmJjY6lTp44CrYiIOJfYWBg6lOujRtEA2HFjs69vZ1aunMbDD2tFBpGsIk1WKMufPz+tW7dOi1OJiIikrQMHoHNnDu7YQTPgzI3N1aoNYd26D/D2tjiwOBFJa/oLGBERyZoMAyZOJC4wkI47dlAZe7B1o127AWzbNkzBViQLSpORWxEREacSEgLdu/PHihW8BBy5sdnPrw5r1vxAtWqFHVmdiKQjhVsREclalizhxIsv8s6VKyxM3OhBYOCTbN48Fw8PdwcWJyLpTW0JIiKSNYSFca1zZwa1bUv5xGBrJW/eV9m48RQ7dixUsBV5AGjkVkREMr2E335j1nPPMfjKFUJubLPSiJde+ZIJEyrj5ubQ8kQkAynciohI5hUZyZrOnem3ZAm7bmzypDA5Ckxg0aLW1K2rB8ZEHjRqSxARkUzph48+olmePDS9EWyz4Y47I+jY/RiHD7dRsBV5QGnkVkREMpW4iAieqVqVpUfMORBcgYK05nquacyYkQdNuy7yYNPIrYiIZBoLRo+mUu7cicE2J9nwZC0PPbGE/fsVbEVEI7ciIpIZxMcz9LHHGP777wD4AsXpwRGvaXz+ObzyCljUhSAiKNyKiIiTSzhwgPeaNmX0uXMA+OKOC2vwqFmPHXOgbFkHFygiTkVtCSIi4pwSErj68ce0CghIDLbleZRwInh9SD02bFCwFZHbaeRWREScz6FDHOrYkad27CAY8MCKN2Mxyr3OX7OhVi1HFygizkojtyIi4jwSEmDMGJZXrkytG8HWm5zEsIWub77Ojh0KtiJybxq5FRER53DoEEb37ny6aRMDAANwpya5/H/ml9n5adTIwfWJSKagkVsREXGsG6O1UVWq0HXTJvpjBlt4iedfXM++fQq2IpJ0GrkVERHHOXQIunfn7KZNtAH+BsAFX9+xzJ3biyee0PxeIpI8GrkVEZGMd2O0lsBAFm7aRDUsN4JtLpo0+Y3jx3sr2IpIimjkVkREMtb+/dCjB7bNm+kBzALAwMWlImPHLqV371IOLU9EMjeFWxERyRixsTBqFHz0ERvj4vgfVvZiA8DdPQ9//x3EQw8VcHCRIpLZqS1BRETS399/Q/XqXBg6lA5xcdSDG8HWiyeeGEBExDkFWxFJExq5FRGR9BMZCUOGcPnzzxljGIwDrt/4yMXFl/nzl/H00/UdWaGIZDEKtyIikj5+/53L3bsz9p9/GAeE39hstVbn+ef7MG1aJ9zcXBxZoYhkQQq3IiKStq5e5cLrr/PFnDlMAK4lfvAQ5csP46efWlO2rGZCEJH0oZ5bERFJM2dnzuStIkUoPmcOIzGDrZUA3N0X8vnnO9i7t42CrYikK43ciohIqm1atoznO3Tg5PXrxN/Y5kNZIhhN3fqtmDnTSunSDi1RRB4QGrkVEZEUi4mKYnKXLjR78kmO3Qi2xSiMGz+SkO0gX37Zmj/+ULAVkYyjkVsREUmRLQsW8FinToQlJADghwu+DOIfhtGwIUyfDqW0HoOIZDCN3IqISPJER7Prf/+jbYcOicG2pqUR14jgktcwJkyANWsUbEXEMTRyKyIiSbd2LYs6daLruXNEArks7rgY0/nbeJ5GjczR2pIlHV2kiDzINHIrIiL3d+kStu7dGda4Mc/cCLYFLFW4bJwj0vt5Jk6EoCAFWxFxPI3ciojI3RkGfPcdkX378sLFi/xwY7M7rxFijKNxY1emT4cSJRxapYhIIqccuZ04cSLFixfH09OT2rVrs2XLliQdN2/ePCwWC23atEnfAkVEHgRHj0Lz5pzp0oVHbgRbC67AdNx9JvLVV66sXq1gKyLOxenC7fz58+nXrx9Dhw5l+/btVKlShebNm3P+/Pl7HnfixAnefvttGjRokEGViohkUTEx8OGHEBDAtlWrqAVsAyAPBmto2vRF9uyBnj3B6nT/FxGRB53T/bb0+eef8/LLL9O9e3cqVqzI5MmT8fLyYsaMGXc9JiEhgc6dOzNs2DBKquFLRCTlVq+GypVhyBAWR0dT32LlLAAV8fH5i6+/bsCqVVC8uGPLFBG5G6cKt7GxsWzbto1mzZolbrNarTRr1oxNmzbd9bjhw4eTL18+evTocd9rxMTEEB4efstLROSBd/YsPPccPPooxuHDjPDy4Wkg2rABj/P44xs5cKAkr7wCFq2eKyJOzKnC7cWLF0lISCB//vy3bM+fPz8hISF3PGb9+vVMnz6dqVOnJukaI0eOxM/PL/Hl7++f6rpFRDKt+HgYNw7Kl4f58wkHynjmZlBkBACenr357rufWb7cjyJFHFuqiEhSOFW4Ta5r167RpUsXpk6dSp48eZJ0zMCBAwkLC0t8nTp1Kp2rFBFxUps3Q82a8OabcO0avxWtQB6LL0ejLwFQvfrnnDo1nk6dXDVaKyKZhlNNBZYnTx5cXFwIDQ29ZXtoaCgFChS4bf+jR49y4sQJWrVqlbjNZrMB4OrqyqFDhyj1nyVyPDw88PDwSIfqRUQyicuXYeBAmDoVDIOdPj60tuXg5MkDN3Zwo0OHd5g3r69DyxQRSQmnGrl1d3enevXqBAUFJW6z2WwEBQVRp06d2/YvX748e/bsYefOnYmvp556isaNG7Nz5061HIiI3Mxmg1mzoFw5EqZMYalh0ChXYapGRHAy8jQApUs/x4kTIcyb97FjaxURSSGnGrkF6NevH926daNGjRrUqlWLsWPHcv36dbp37w5A165dKVy4MCNHjsTT05OAgIBbjs+RIwfAbdtFRB5o27dD796c2LSJmcAMV1dOx8fD5TOABVfXEnTv3oPJkwdg1fxeIpKJOV247dChAxcuXGDIkCGEhIQQGBjIihUrEh8yO3nypH7jFRFJqkuXuPbOOyyaOZPZwO/27fHxQG4slhf43/9e5YsvSuHp6bgyRUTSisUwDMPRRThSeHg4fn5+hIWF4evr6+hyRETSROS1a7zbqhVHN2zgj/h4ohI/sQBNgB5UqdKOWbM8CAx0VJUiIkmTnLzmdCO3IiKScjabja3Tp9OuVy/OxMUlbi+Yw5+L118hLq4Lnp7F+OADeOstcNX/BUQki9FvayIiWcTJrVtp1LgxxyPMOWotQK0ipbjm/Q37D9UFLDzyiDlJQtmyDi1VRCTdqHlVRCSzi4vjwrBhtKhdOzHYti5ajDe77mLruSPsP1QPPz8LkyfD778r2IpI1qaRWxGRzGzVKoJ79qTFsWMcAzwtFl5t+xY/bv+Un2abuzz7LIwdCwULOrJQEZGMoXArIpIZHT4Mb73Fhp9/5ingMlAsVx4CHv6DLxZXBMDfHyZNgiefdGilIiIZSm0JIiKZydWr5pNglSqx+OefaYoZbEsVr8HVhH0sW14RqxX69oX9+xVsReTBo5FbEZHMICEBpk2DwYPh4kUmAG8ABpA7d2uOnpgLeFG1qvnAWPXqji1XRMRRNHIrIuLs1qyBqlWhZ0+Mixd5L3duXscMtlbrq1y6tAgvLy8++wy2bFGwFZEHm0ZuRUSc1ZEj8Pbb8NNPAMTnyMH/ypdnxubNN3b4EJttEI8/bmHSJChRwnGliog4C4VbERFnExYGH30E48ZBXBy4uBD1yis8ffg4v65egfmXblPIl68H48ZBhw5gsTi6aBER56C2BBERZxEfD5MnQ5ky8NlnZrBt3pwrf66n+tp9N4KtB/AjPXr04MABeO45BVsRkZtp5FZExNEMA375Bd59Fw4eNLeVKweff87SuHw8/UhT4uPDAV/8/X/m228foWFDh1YsIuK0NHIrIuJI27ZBkybw1FNmsM2dG778koiN23hs7C5at6lzI9ha6dFjNcHBCrYiIveicCsi4gj//APPPw81asDateDhQfTbbxM0YwaPL9uGb548rFr1HhCPm1s+pkz5hWnTauLp6ejCRUScm9oSREQy0tWrMHIkjBvHpZgY/gY2Va7Mnz4+bBw/npjPPrtl9xdfnMzUqa9gtaqxVkQkKRRuRUQywNGDBxnxwgt47NxJaEwMO4Dj9g/37Llpz4JAMypXzsmMGb2oUaNshtcqIpKZKdyKiKSjq5cv89uHH9Lryy+5aLPd9rm/f2mioupw8WI9oCFVq5bj668t1KyZ8bWKiGQFCrciIunBMFg5ciRtBw8myjAA8zfcGiVK8Myrr1KuQnWWLavKlCk5sdkge3b4+GN47TVwcXFs6SIimZnCrYhIWtuyhZN9+/Lyxo1EARagb506fLBwIT4FC7FoEfTsCWfOmLt36ACffw6FCjmyaBGRrEHhVkQkrRw8CIMHc2zRIpoAp4A8np4s+/FHaj3+OMeOQYcn4Ndfzd1LlYKJE6F5c0cWLSKStWgqMBGR1Dp9Gl56CSpVInjRIh4B/gHKlizJ9uBgqjR+nI8/hkqVzGDr7g7vv28+R6ZgKyKStjRyKyKSUpcvm9N6jR8PMTEcAhp7eHAuJoaKFSsSFBTEgQMFePRROHTIPKRJE5g0yVyATERE0p7CrYhIckVEwJdfwujREBYGwKHq1Wl04gQhly4REBDA998H8e67+fj2W/OQfPngiy+gY0ewaMpaEZF0o3ArIpJUkZHmsOsnn8DFi+a2hx7icO/eNB469EawrUzHjkHUr5+XsDAzyL76qjkTQo4cDq1eROSBoHArInI/0dEwZYrZghASYm4rXRo++IBjtWvTuFEjzp07R8mSlbHZ1jBoUB4Aqlc3s3CtWg6sXUTkAaNwKyJyN7GxMGOGOex6+rS5rXhx82mwrl05efYsTR55hDNnzuDnV4Fjx1YDeciZE0aMgJdf1py1IiIZTeFWROS/4uNh9mz48EM4ccLcVqQIDB4M3buDuzshISE0bdqUf/75B6u1DGFhQUA+XnrJHODNk8eRX0BE5MGlcCsiYpeQAN9/D8OGwZEj5rYCBeC998xhWE9PAC5fvky9eo9y7NgRoDg2WxDVqxdk4kSoXdtx5YuIiMKtiIgZan/4AYYPhwMHzG158sCAAebTYF5eibtu3XqCBg1qER19ASiIr28Qn3zirxYEEREnoXArIg+shJgY+P57XEaN+nci2pw54d13oXdv8PFJ3DcmJoEXXpjBvHmvAzGAC+3a/cbXX5dUC4KIiBNRuBWRB058ZCT9WrVi6u+/U84wqAL4eHiQrXp1XGvVgqtXif/gAyIjIwkPD2f37hPs3bsDw4gEwGLx4t13P2fUqADHfhEREbmNwq2IPDBir13jmXr12LV/PycTEgDYdeNFTAxs3Gi+7srCM898zuzZr5Etm3sGVCwiIsmlcCsiWV9UFEydSvv+/fk5OhoAX+CpatWo2bEj0TYb165dIzo6moSEBGw22LvXlfXrsxETkx3ISYUK2xk6tBUdOjzu0K8iIiL3pnArIlnX9esweTJ8+ikzQ0NZemNz68qVmRMUhE/evLcd8uef8PrrsGuX+b5aNZg4ER5+uEfG1S0iIilmdXQBIiJp7vJlc47a4sXh7bdZERrKyzc+6v/22yzZvfu2YHv6NHTqBI88YgbbnDnhq69gyxZ4+OEM/wYiIpJCGrkVkazj1Cn44gtzqdzr1wHYWaQI7S9cICEmhi5dujBy9OhbDomOhs8/Nxchi4wEiwVeeQU++kgLMYiIZEYauRWRzO/AAXPlsJIlzXB7/TpUqcKZiRN5wmYjIiaGJk2aMG3aNCwWCwCGAUuXQqVKMGiQGWzr1YOtW81OBgVbEZHMSSO3IpJ5bdoEn3wCP/3077ZGjWDAACLq1qVVw4acPXuWihUrsmjRItzdzRkODh6EN9+ElSvNQwoVgk8/hY4dzZFbERHJvBRuRSRzMQz49Vcz1K5bZ26zWKBtW+jfH2rVwmaz0eXpp9mxYwf58uVj2bJl5MiRg/BwcxGyceMgPh7c3eGtt8zVdW9ar0FERDIxhVsRyRxiY2HBAnOIdfduc5ubG3TpAu+8A+XLJ+76/vvvs2TJEtzd3VmyZAlFixZn1ixzNd3QUHOfVq3MXtvSpTP+q4iISPpRuBUR53b5Mnz9NUyYAGfPmtt8fOB//4O+faFw4Vt2nzdvHiNGjABg+vTpuLjUoU4dc9YDgLJlYexYaNEiA7+DiIhkGIVbEXFOwcFm/8CsWebTXgAFC0Lv3vDqq+ZcXf+xY8cOXnzxRQB69XqXNWuep0sX8zMfHxgyBPr0MdsRREQka1K4FRGnYBiG2U/7xx9YvvgCfvnFfA8QGAj9+kGHDndNphcuXKBNmzZERUVRvnwLZs8ewbVr5mddu8KoUWY2FhGRrE3hVkQcLjYigpFduvDJTz/RxTD4HPAGszG2b19zBoR7TGMQHx9Px44dOXnyJG5upTl4cC7gQvXqMH481KmTMd9DREQcT+FWRBzm3K5dtHz8cU6HhnLxxijtFGCumxsdWrem82uv8cgjj+Byn/m5XnjhDYKCggBv4uKWkDdvDkaONKe+tWo2bxGRB4p+2xeRjGUY5vy0zz9Pj8BAdoaEcNEwyAM0L1OGEkWLEhEXx/SFC2nSpAmFChWiR48eLFq0iEuXLt1yqnPnIihfvj3fffcVABbLDN58sxLBwdCjh4KtiMiDyGIY9qa2B1N4eDh+fn6EhYXh6+vr6HJEsq6oKJg3z5z1YPt2fgeaAgbwWtOmfP7jj3hkz45hGPz55598++23LFq0iCtXrtxymkqVKlGzZi2CgxPYtGkZhmEG3mzZSrF16xEqVszwbyYiIuksOXlN4VbhViR9HT8OX30F06eb03oBF93dqeLqytnISHr06MG0adPueGhcXBx//PEHv/zyC7/99hsHDhy4bR9X1xI8/XQXvvlmEB4emgZBRCQrUrhNBoVbkXRgs8Hq1eYo7c2zHhQrhvHqq7RZu5alK1ZQvnx5tm7dire3931P+c8/0KfPeX76aROwHpiEv38Ae/b8gZ+fZ7p+HRERcazk5DV1pIlI2rlwAcaMgQoVoHlz+PlnM9g+9hgsXQpHjzIlRw6WrliBu7s733///X2DbUQEvP++uQDZTz/lw2Jpzcsvf8rp02GcPPmXgq2IiNxCsyWISOrYbPD77zBlCvz4I8TFmdt9fc3pCl57zVwWDAgODqZv374AjBo1isDAwHueds4cGDjw34XJGjY0VxczD9NvXyIicjv930FEUiYkxFw9bOpUOHbs3+01asDLL0OnTuayYDfEx8fTpUsXoqKiaNasGX369LnrqTdtgjff/HfJ3BIl4LPPoG3be053KyIionArIsmQkACrVpmjtD//DPHx5nZfX3j+eTPU3mU0dtSoUWzZsoUcOXIwc+ZMrHeYp+vUKRgwAObONd/7+MCgQWbQ9VT3gYiIJIHCrYjc3+nTMGOGOePByZP/bq9TB155Bdq3h3v0zu7cuZNhw4YBMH78eIoUKXLL55GRMHq0+YqKMkdnu3eHjz+GAgXS5RuJiEgWpXArIncUfu4cvdq3532bjbKbN/8740HOnNC1K7z0EgQE3Pc8sbGxvPDCC8THx9OuXTs6d+6c+JlhmKO0AwaY+Rmgfn2zr7Z69XT4UiIikuUp3IpIIlt8PENffJGQv/5idnAwsZiTbu0E/Bo2NNsOnn46WT0CI0aMYNeuXeTOnZuvvvoKy42m2S1boE8f2LzZ3K9YMfj0U3jmGfXViohIyinciggcOQLffsuLY8bwzfXriZstwAmgbpky/DRtGqVLl07WaXfv3s3HH38MwIQJE8iXLx9nzpgzIHz7rbmPt7f5vl8/yJYtbb6OiIg8uDTPrciDKizMnOmgfn0oU4ao4cMJuhFsK+XIwW9jxrBlyxYKFSrE/sOHqVmzJsuXL0/y6ePj4+nRowfx8fG0adOGJ5/swAcfmLOC2YNtt24QHGw+NKZgKyIiaUEjtyIPkqgoWL4c5s0zVw6Ljja3W618XrIkp48coVDBgmwODsbnxjReW7dupW3btvz111888cQTDBgwgOHDh+Pm5nbPS3355Zds3boVPz8/GjacRLlylsT5auvWNftqa9ZMx+8qIiIPJC2/q+V3JauLizOn75o3D5YsgWvX/v2sYkXo1o0zTZpQtmFDIiMj+e677+jUqdMtp4iJiaFfv35MmjQJgFq1avHtt99S9sbiDP914sQJKlWqRGRkJP7+Uzl16iUAihc3Z0RQX62IiCSHlt8VedAlJJirhr3yijmX1hNPmL0A166Bvz+88w5s2wZ798K77zJg3DgiIyOpV68eHTt2vO10Hh4eTJw4kQULFpAjRw62bNlClSpVGDNmDPH2uW5vMAyD5557kcjISKAhp071wNfXDLUHDpizhinYiohIetHIrUZuJaswDPjrL/j+e1iwwFxBzC5/fnj2WXjuOXj4YbhpAYUtW7ZQu3ZtAP7++29q1Khxz8ucOnWK7t27ExQUBEDVqlWZMGECdevW5cSJcJo27cqxYz8Brlite3j11fIMHQp586b5NxYRkQdEcvKawq3CrWRmCQmwYQP8+CMsXnzrAgs5c5rTdj33HDRsCK63t9gbhkGDBg3YsGED3bp1Y9asWUm6rGEYzJgxg3feeYcrV64AUKRIdU6fPg5cBsDHJ4AtW/ZQoUJqv6SIiDzoFG6TQeFWMp3YWFizxgyzS5bAhQv/fubjA61bQ8eO8Oij4O5+z1MtXryYp59+mmzZsnH48GEKFy6crFJCQ8/TufN7BAXNAMzfStzcStOxYxcmT36HbJoCQURE0kBy8ppmSxBxcoZhEHXxIjMHDuTjuXMZbxg8bZ/lAMwR2latoF07eOyxJM+pFRcXx4ABAwB46623kh1st22Dfv3ysW7dNKAL0JLAwEfZtOkHPD3vPZOCiIhIelG4FXFWZ8/y/Qcf8MqMGcQkJBB3Y3MfoG2BAljbtjUDbcOGcJ9pue5k6tSpHD58mLx58/LOO+8k+bjTp815aWfPNt97esJbbzXknXeukT07WK16TlVERBxH4VbEWdhs8PffsGyZ+dq+nVlAxI2PC7m4cMEwOGOz8cVbb/HW22+n+FIREREMHz4cgCFDhiSpJSciwlwe99NPzelyAZ5/HkaMMCdg0OQrIiLiDPR/IxFHCguDH36AF16AggXNmQw+/BC2b+cqsOHGKOiXAwZwOjaWCV99BcDA995jx44dKb7s2LFjCQ0NpVSpUrzyyiv33DchAWbMMFcWGz7cDLb168OWLebsYmawFRERcQ56oEwPlElGMgzYtw9WrjRHZ//8E26eJ9bXF5o3hyefZFRwMAM//piAgAB2796NxWLBMAzatWvHkiVLKFu2LFu3biV79uzJKuHSpUuULFmS8PBw5s6de8d5be3WrIF+/WDXLvN9yZLmfLXt2mmuWhERyTh6oEzEmZw5A6tX//u6ef5ZgPLlzUUWnnjCHBJ1cyM2NpYvixcH4O2338ZyI0laLBamTZvG33//TXBwMD179mTOnDmJnyfF6NGjCQ8Pp0qVKnTo0OGO+xw6ZK7z8PPP5ns/PxgyBHr1Ag+PZN8BERGRDKNwK5LWwsPhjz/MJW9XrzaX5bpZtmzwyCPQsqUZaEuVuu0U33//PefOnaNQoUK3jazmzp2b+fPn07BhQ+bOnUu9evV47bXXklRaSEgI48ePB+Cjjz667eGv8+dh2DD4+muzHcHFBV57zQy2efIk4x6IiIg4iMKtSCoYhoElJsZ8ECwoyAyzmzebydDOaoUaNaBZM3Pu2Tp17jn8aRgGY8aMAaBPnz6432Gu2nr16vHJJ5/w9ttv8+abbxIYGEjdunXvW+/IkSOJioqiTp06PPHEE4nbIyPhiy/gk0/MFXrBnF1s9GhzYFlERCTTMJzQhAkTjGLFihkeHh5GrVq1jL/++uuu+06ZMsWoX7++kSNHDiNHjhxG06ZN77n/f4WFhRmAERYWlhalSxYXGxtrBO/YYfzx6adG+woVDB+r1fjK1dUwzG7af1+lSxtGz56GsWiRYVy+nKxrrFq1ygAMb29v48qVK3fdz2azGe3btzcAo0CBAsbp06fved5Tp04Z7u7uBmCsXr3aMAzDiI83jJkzDaNw4X9Lr17dMH7/PVkli4iIpKvk5DWnG7mdP38+/fr1Y/LkydSuXZuxY8fSvHlzDh06RL58+W7bf+3atXTs2JG6devi6enJJ598wmOPPca+ffuSPSm9yB1duADr18Off/La7NlMu3Tplo/H22z0zJcPGjUyR2abNYMb/bIpMXbsWAC6d+9Ojhw57rqfxWJhxowZHDx4kD179vDUU0+xbt06vL2977j/iBEjiI2NpWHDhjRp0oTffjP7anfvNj8vVsyc1uu558zBZhERkczI6WZLqF27NjVr1mTChAkA2Gw2/P39ef311xNXU7qXhIQEcubMyYQJE+jatet999dsCXILm818muqvv2DTJnM2g5t6ZusBG4FsFgv18+dnVUgILlYrR44coXiJEqm+fHBwMOXKlcNisXDo0CHKlClz32OOHz9OrVq1uHjxIm3atGHhwoW4uLjcss+pU6coVaoUcXFxTJv2OwsWNOK338zPcuQwF2Xo3dtckEFERMTZJCevOdX4TGxsLNu2baNZs2aJ26xWK82aNWPTpk1JOkdkZCRxcXHkypUrvcqUrOTsWViyBAYOhKZNzaRXsSJ07w5TpvwbbCtV4tqLL7LjRv/r2s2b+e3cOZo1a0aCzcbnX3yRJuXY/1DXsmXLJAVbgBIlSrBkyRI8PDxYsmQJffr04b9/Zh01ahRxcXEUKNCIl182g62bG/TtC0eOwNtvK9iKiEjW4FTh9uLFiyQkJJA/f/5btufPn5+Q/06fdBf9+/enUKFCtwTkm8XExBAeHn7LSx4Q4eHw++/mU1Pt2kGRIlC4MLRtC6NGmZO6XrsGXl7QoIGZ+JYsgYsXYe9eljRqRFRsLGXLlqVmzZqA+esNYNq0aZw/fz6V5YUzc+ZMwHyQLDnq1avH7NmzsVgsTJw4kY8++ijxsx07jvDVV18DEBIyFMOAZ581c/vnn0Pu3KkqW0RExKk4Xc9taowaNYp58+axdu1aPO8yDDVy5EiGDRuWwZVJhjIMOHnSXHlg585//3ns2O37Wq0QEAC1akHt2uY/K1YE19v/0/juu+8A6Ny5c+K8sk2bNqVWrVps2bKFL774gpEjR6a47G+++YaIiAjKly9/1z+c3cuzzz5LaGgob7zxBkOGDCFbNh927crJnDm9gQSgEPXqNWTMGPOrioiIZEVO1XMbGxuLl5cXCxcupE2bNonbu3XrxtWrV/npp5/ueuxnn33GRx99xOrVq6lRo8Zd94uJiSEmJibxfXh4OP7+/uq5zaROHjrEitmzcT17lu+Dgvj77FlGubnRMzr6zgcULfpviK1dG6pVg7s8gHWz0NBQChUqhM1m4/Dhw5QuXTrxs6VLl9K6dWuyZ8/OiRMnUtQSY7PZqFChAsHBwUycODHJ89beydChHzB8+H//AOfKM8/0Z8GCj7SymIiIZDqZdoUyd3d3qlevTlBQUGK4tdlsBAUF0bt377seN3r0aD7++GNWrlx5z2AL4OHhgYeWWMp8Ll82/x79P6+XT5zgt//s+n1CAj1dXc0R2MBAqFLl33+m8O/gFyxYgM1mo1atWrcEW4BWrVpRpUoVdu3axbhx41L0NwNBQUEEBweTPXt2unTpkqIawXwObs2aocAW4FfASqtWH/HNN6+TI4e3gq2IiGR5ThVuAfr160e3bt2oUaMGtWrVYuzYsVy/fp3u3bsD0LVrVwoXLpz417+ffPIJQ4YMYe7cuRQvXjyxN9fHxwcfHx+HfQ9JgZgYOHHCbB84fPjWIHuXflb7pFzFPTwoX6AAK/75hxMFCmAcP44lDZ+Qmjt3LgCdOnW67TOLxcL777/PM888w7hx43jzzTfJmTNnss4/adIkwPz1nT179mTXd+wYvPcezJ8PYMHdfSn+/q35/vvPqFmzQrLPJyIiklk5VVuC3YQJE/j0008JCQkhMDCQL7/8kto3mgQbNWpE8eLFmTVrFgDFixfnn3/+ue0cQ4cO5YMPPrjvtTQVWAYyDPPhrGPHzNfRo//+fOwYnD5t7nM3RYpAhQqJr7gyZcjdpg3XIiL4+++/qVChAnnz5iUqKort27dTtWrVNCn7+PHjlCxZEqvVyunTpylYsOBt+9hsNgIDA9mzZw+DBw/mww8/TPL5T58+TbFixbDZbOzbt4+KFSsm+djLl+Gjj2DCBIiLA4sFunWDDz80b5eIiEhWkGnbEux69+591zaEtWvX3vL+xIkT6V+Q3J9hwNWr/L1iBWtXr8YvNpY1O3aw8/RpPi5RgqcTEsxRWfvarnfj7Q2lSkHJkua6r/YwW748/GdEc9O6dVyLiCBv3rxUq1YNq9VKy5YtWbRoET/88EOahdsFCxYA5h+s7hRswZyy7oMPPuDpp59m7Nix9OnThzx58iTp/FOnTsVms9GwYcMkB9voaJg40Qy2V6+a2x59FD791Oy+EBEReVA5ZbgVJ2EYcP262RJw4cK/L/v70FA4c8YccT1zBiIjeRtY95/TLNy5k6dv3lC48L8B1v6yv8+bl6Q2hq5YsQKA5s2bY72xpFb79u1ZtGgRCxYs4OOPP06c1SA15s2bB8Bzzz13z/3atm1L1apV2bFjB5988gmffvrpfc8dFxfH1KlTAXj11Vfvu7/NZrYevPee+WcFgMqVzVDbvPl9DxcREcnynLItISNl6bYEw+DCmTNEXLiAW1wcJ4OD2b5tGy0rVqSkl5c55Hf1KoSF/fvz1avm33XbA+zdZh24i0AXF3YlJFAqWzayZ8vGzsuXaVK+PEFffGEuSVu8eJqtFlCtWjV27NjBt99+y/PPPw9AREQE+fLlIyoqiq1bt1K9evVUXePQoUOUL18eV1dXQkJCyH2fB9KWL1/OE088gaenJ0eOHLnvEtBLliyhbdu25MuXj1OnTuF+Y5GIO/njD3Pq3a1bzfeFCpkjt127wn8WJBMREclSMn1bQla2cfJk9m/cSEJ8PLtPneL0lSu817gxtQsXhvj45L3i4iAqCiIj/31dv37L+9Y2G/9d2+0U8Elyivb0hHz5zFHVvHn//TlfPnMU9sZiCEahQoSULAmhocxevRrDMKhfvz57Ll3CaN48TUZR7UJDQ9mxYwcAjz32WOJ2Hx8fnnzySX744QfmzZuX6nD7ww8/ANCsWbP7BluAFi1a0KBBA/7880+GDRvGlClT7rn/11+biyt07979rsH2wAEYMACWLjXf+/iY7/v2NdebEBERkX9p5DaDR25fCwjgq337btk2GEj640fJUwnYj7kUnSsQCzyXMyffV6tmLjVrf/n5/ftzzpy3hllv7yS1Chw7doxSpUrh5uZGeHg4VquVHDlyEBUVlewHpe5nzpw5dOnShapVq7J9+/ZbPlu0aBHPPPMMRYsW5fjx44ktCynx0EMPsWfPHmbOnMkLL7yQpGM2bNhA/fr1sVqt7N27lwoV7jxbwYkTJyhZsiSGYXDkyBFKlSp1y+dnz8IHH8D06WY7gosLvPIKDB0K/1nET0REJEvTyK0Te+ihh2h18SJWq5VtV69yOiqK8PLloW5dc1Us+8vF5db3d3t5ed368va+5b3RqBEEB/PL8uUcPHiQfv36EfXII+aysmls0yZzjLhatWqJK8TVqVOHNWvWsHbt2jQNtytXrgTMftv/atmyJdmzZ+fkyZNs2rSJevXqpegaBw4cYM+ePbi5udG6deskH1evXj3atGnDkiVL6N+/P0vtQ67/MWPGDAzDoFmzZrcE2/BwGD3aXBo3Ksrc1rq1uUJw+fIp+ioiIiIPDuMBFxYWZgBGWFhYhl976NChBmD06NEjXc4fFRVluLi4GIBx5swZY+3atQZgFC1aNF2u16tXLwMw+vbtm7ht+PDhBmC0b98+za6TkJBg5M+f3wCMNWvW3HGfrl27GoDx2muvpfg6w4YNMwCjZcuWyT72wIEDiff+TjXGxcUZhQsXNgBj/vz5hmEYRkyMYYwbZxh58hiG+TSfYdSpYxjr16f4K4iIiGQJyclrKf/7Wkm18jeG4Q4ePJgu59+/fz8JCQnkzp2bggULEhgYCMDJkye5dOnSvQ9OAfvIbZ06dRK3NWrUCIA//vgDI406YPbs2UNoaCheXl7UrVv3jvt07NgRMHtm4+PjU3SdhQsXAuYMDMlVvnz5xNkP+vbtS0JCwi2fr1y5kjNnzpA7d25atWrN99+bo7J9+phTAZcrB4sXw4YNkMKBZxERkQeSwq0D2cPtoUOH0uX8u3btAqBKlSpYLBb8/PwoWbIkADt37kzTa12/fj3xejeH21q1apEtWzbOnz/P/v370+Raq1atAqBhw4Z3XUq5adOm5M2blwsXLiTunxyHDh1iz549uLq68tRTT6Wozg8++IAcOXKwa9cupk2bdstn06dPB6Bhwy40aOBBp05w/DgUKABffw1790LbtkmeFU1ERERuULh1oDJlygBw8eJFLl68mObn3717N2D2+drZFzb470NYqbVt2zYSEhIoXLgwRW5aGsvDwyOx5/X3339Pk2vZw+rNsyT8l5ubGx06dADgu+++S/Y1Fi1aBJghOVeuXCmoEnLnzs3w4cMBGDRoEJcvXwbg/PnzLF36MwCLF/dg2zZzfYoPP4QjR8yHxlzVDS8iIpIiCrcO5O3tTdGiRYH0Gb21j6TeHG6rVasGpH243bx5MwAPP/zwbZ81btwYSJtwGx0dzbp15jIRjz766D337dKlCwA//vgjERERybqOPdw+/fTT99nz3l599VUCAgK4dOkSgwYN4u+/z1GpUmsSEuKB0ri5BfD66+ZKxIMHm88DioiISMop3DpYevXdGoZxx5Hb9Aq3d+q3tbOH27Vr12Kz2VJ1nY0bNxIdHU3BggXvO/tCzZo1KVOmDJGRkfz4449JvsaJEyfYvn07VquVNm3apKpeV1dXJkyYAMDkyZOpVasEFy+afxDIl68oBw7Al1+aM66JiIhI6incOpg93B44cCBNzxsSEsKlS5ewWq23hEB7uA0ODiY8PDxNrmUYxj1HbmvUqIGPjw+XL19ODNwptXr1asActb3fohAWiyVx5bLZs2cn+Rr2IPzII4+QN5WpMyoKNm9uiNVqn+s2Bh+fWrz33tcEBy/mP1PbioiISCop3DpYeo3c2kNk2bJlyZYtW+L2fPnyJfbEptVDZSdPniQkJARXV9fE8HwzNzc3GjZsCEBQUFCqrmUPt02bNk3S/vbWhKCgIE6fPp2kYxYvXgxAu3btUlChKSEBZs6EsmXN1cRstnlYLNlo334AYWGb+fjjV/Dz80vx+UVEROTOFG4dzL56VVqH2z179gBQuXLl2z6zL0mbVq0Jf/31F2DOynBzkL5ZkyZNAFizZk2Kr3PlyhW2bt0KJD3clihRgkceeQTDMJgzZ8599w8NDWXDhg0AKWpJMAxYtgwCA+HFF+H0afD3h2++eYiYmAgWLBiJ1aopEERERNKLwq2D2Udujx8/TnR0dJqd9079tnZp3XdrD7e1a9e+6z72cLtu3Tri4uJSdJ3ff/8dwzCoUKEChQsXTvJx3bp1A2DmzJn3nWt36dKlGIZBjRo18Pf3T1Z9W7ZA48bw5JPmVF45c8Knn0JwMHTtCm5u+s9NREQkven/tg6WP39+cuTIgc1mIzg4OM3Oe69wax+53bZtW5pcKynh9qGHHiJ37txERESwZcuWFF3H3tKQ1FFbu/bt2+Pl5UVwcHDig293Y++3bdu2bZLPf+QIPPss1K4Nf/wBHh7wzjvmDAhvvw03ViIWERGRDKBw62AWiyXN+27j4uISH1C7V1vCwYMHkz1F1p2uZQ/J9wq3Vqs1MZTa+2aTy97SkNxwmz179sRVxmbMmHHX/cLDwxMDdFJaEkJDoXdvqFABfvjBXHChWzdzpHb0aHPkVkRERDKWwq0TsPfdptWMCcHBwcTGxuLj40OxYsVu+7xAgQIUKlQIm82W6ofK9uzZQ3R0NDly5EhclOJumjVrBqQs3J49e5aDBw9isVgSH05LjhdffBGA+fPnc+3atTvus3LlSmJjYylTpkziv5M7iYiAYcOgdGmYOBHi46FFC9i5E2bNghtTF4uIiIgDKNw6gbQeubU/TBYQEIDVeud/xTVq1ABS35pgb0moWbPmXa9lZw+3mzdvvmvAvBv7AhDVqlUjZwqGRBs0aEDp0qWJiIhgwYIFd9znp59+AqB169Z3nGYsLg6++soMtR98YIbcmjVhzRpYvhzu0AEiIiIiGUzh1gmk9cjtvWZKsLO3JthnH0gpe//svVoS7EqUKEGpUqWIj4/njz/+SNZ17C0J9gfTkstisfDSSy8BMG3atNs+j4uLY9myZYAZbm9mGLBoEQQEwGuvme0IpUrB/Pnw11/mQ2QiIiLiHBRunYA93B46dIiEhIRUny854Ta1I7f2cFurVq0k7W8fvV21alWyrmMPt41TkSS7deuGq6srmzdvTrxHduvXr+fq1avkyZPnllXW/vwT6taFZ54xe2nz5oXx42H/fvMhsvusIyEiIiIZTOHWCZQoUQIPDw+io6M5ceJEqs+XnHB78ODBZLcI2IWHhyeONtesWTNJxzz66KNA8sLtiRMnOHHiBC4uLtSvXz/5hd5QoEABnnrqKQC+/vrrWz5bunQpAE8++SQuLi7s2wdPPQWPPAKbN4OXF7z/vjkzQu/e4O6e4jJEREQkHSncOgEXFxfKli0LpL414dq1a4kB+V7htkCBAhQpUgTDMNixY0eKrrVt2zYMw6Bo0aIUKFAgScc0adIEq9XKgQMHOHXqVJKOWbt2LWAG6OzZs6eoVruePXsC8O2333L9+nXAXD74559/BuDhh1vRvbvZP/vzz+DiAv/7nxlqhw8HX99UXV5ERETSmcKtk6hYsSKQ+nC7b98+AAoWLEju3Lnvua/9obKU9t0mtyUBIGfOnIn7J3X01v4wWWpaEuyaNm1KyZIlCQ8P5/vvvwfM0eujR4/i4uLOG288xqxZYLNBu3awbx9MngwFC6b60iIiIpIBFG6dRFo9VLZ3717AnCnhfuytBH///XeKrmU/LjnhFv5tTfjtt9/uu69hGIkjt2kRbq1Wa+Lo7aRJk7hwIZLWrd8CICGhGrGxPjRubD4otmgRlCuX6kuKiIhIBlK4dRL2cLt///5Unccebu/VkmCX2pFbe7hNar+tXfPmzQFz5PZ+D9CdOHGCkydP4ubmRt26dVNU5391794dDw8PduzYQf78xTl8+FcAvL1dWbECgoIgmXldREREnITCrZO4uS3BMIwUnyc5I7f2cHvkyBGuXLmSrOuEhoZy8uRJLBZL4sNpSVW7dm38/Py4fPnyfYP1zf223t7eybrOndhssGpVHmw2c3ELw7iAi0sRmjbtyf7939O8uWZAEBERycwUbp1EmTJlcHFxITw8nLNnz6b4PMkJt7ly5aJkyZJA8qcEs4/ali9fPtkPebm6uiZOCbZy5cp77msPt40aNUrWNf7LMGDlSqhRAzp1gri4IYCVgIAnuHLlCKtXf0XRokVSdQ0RERFxPIVbJ+Hh4UGpUqWAlPfdXrhwgdDQUODfkeD7SWnfbUpbEuwef/xxAH799dd77pcW4XbLFmjaFB5/HHbsgOzZ4cMPOxMSEs6ePb+QPbtHis8tIiIizkXh1onYA2lK+27tMyWULFkyyX+Fb38YzD7zQVLZ2wlSG263bNnC5cuX77iPvd/W1dU1Rf22hw6Ziy/Urg2//27OTdu3Lxw7BoMHQ/78qW9zEBEREeeicOtEKlWqBKQ83CanJcHOHm6TM3JrGEaqw22RIkUICAjAZrPdddYE+xK9NWrUSFa/7Zkz8MorUKmSOeOBxQLdupkrjH3+OeTJk6KSRUREJBNQuHUiqR25tYdbe0hOiqpVq2K1Wjlz5kySe31Pnz7N+fPncXV15aGHHkpRrQAtWrQAYPny5Xf8fN26dQA0bNgwSee7cgX694fSpWHqVEhIMFcZ270bZs2CYsVSXKqIiIhkEgq3TsQebvft25eiGRNSMnLr7e2duH9SR2/to7YBAQFky5YtmVX+q2XLlgCsWLECm8122+f2cPvII4/c8zyRkfDJJ1CyJIweDdHRUL8+rF8PP/0EybgdIiIikskp3DqRcuXKYbFYuHz5MhcuXEjWsYZhJPbcJmfkFv5tTfjrr7+StL893NqnEkupevXqkT17di5cuHDbbA1nz57lyJEjWK1W6tWrd8fj4+PNEdoyZWDAALh61QyyP/8M69bBXQ4TERGRLEzh1olky5YtcWoue1BNqnPnznH16lVcXFwol8xltZL7UFlahVs3Nzcee+wxAJYtW3bLZ3/++ScAgYGB+Pn53fKZYcDChWZP7SuvwNmzZsvB7Nmwcyc8+aTmqhUREXlQKdw6mZT23dpbEkqXLo2np2eyjr35obI7tQfc7OaHyVIbbgGeeOIJAH755ZdbtttbEho0aHDL9jVrzNkP2rc3HxDLkwfGjjVnRujSBVxcUl2SiIiIZGIKt07G3lKQ3JHblLYk2I/x8vIiPDycgwcP3nPfEydOcPnyZdzc3JLV23s39r7bbdu2ce7cucTt9pFbe7jdsQOaNzfnq/37b/D2hiFD4OhR6NMHPDRVrYiIiKBw63QcEW5dXV0TR2Hv15pg742tXLkyHmmQKPPnz584nZh91oQrV64kjkQXLlyfjh2hWjX47Tdwc4Pevc1QO2wY+PqmugQRERHJQhRunUxKZ0xITbiFpD9UZg+31atXT9F17uTJJ58E4OeffwZgw4YNGIZBjhxladAgP/Pmmft16gQHD8L48ZA/f5pdXkRERLIQhVsnU758eSwWC5cuXeL8+fNJOsYwjMQe3ZSG29q1awOwefPme+5nD7dp0W9r16pVKwBWrVpFcHAoPXuOAODq1UrEx0OLFmZbwnffmdN9iYiIiNyNwq2T8fLySpwxIakPlZ05c4bw8HBcXV0pW7Zsiq778MMPA7Bnzx6uX79+x30Mw0iXkdvAwEAKFSpMZGQk5cqV5syZTQD4+UWxdi0sXw6BgWl2OREREcnCFG6dUHL7bu37lSlTBnd39xRds0iRIhQuXJiEhITb5py1++eff9L0YTKAuDiYMsVCaKh9ed0I3NxK07btuxw4MJMkLk4mIiIiAijcOiV73639oar7sY/w2o9LKfvo7d1aE7Zv3w6YK5Ol9mEymw3mzoUKFaBnT0hI6A24Ubt2R8LDD7B48ScULFggVdcQERGRB4/CrROyj4omd+Q2pf22dkkNt6lpSTAMcwWxwEDo3Nmc9SBfPvjyy9eJiIhk8+a5eHq6pvj8IiIi8mBTinBCN4dbwzCw3Ge5rbQeud20adMdr2sPt9WqVUvR+deuhffeg01mSy1+fvDuu/DGG+DjA/rlKCIiIqmlkVsnVK5cOaxWK1euXLllYYM7SYuZEuyqVauGq6srISEhnDx58rbr2Htxkxtut241F2Bo3NgMttmywYABcOyYGXbNYCsiIiKSegq3TsjT05MyZcoA929NOHv2LGFhYbi4uCQek1JeXl5UqVIFMEdv/3ud8+fP4+LiwkMPPZSk8x04AM88AzVrmgswuLrCa6+ZrQgjR0KuXKkqV0REROQ2CrdOyj4Ke7+Hyg4cOABA6dKl02TFsLp16wK3h1t7S0LFihXJli3bPc/xzz/QvTsEBMCiRWCxQJcucOgQTJwIBQumukwRERGRO1K4dVL2vtv7hdu06re1q1OnDnD3cFu1atW7HhsaavbPlikDs2aZMyK0aQN79sDs2VqAQURERNKfnuBxUkkduU2vcLtjxw6ioqISR2l37NgB3Lnf9upV+PRTGDsWIiPNbU2bwogRcGNVXxEREZEMoZFbJ3XzjAk2m+2u+6V1uC1WrBgFChQgPj6erVu3Jm63h9ubR24jI2HUKChRwgyykZFmmF292nwp2IqIiEhGU7h1UvbVxq5fv84///xz1/3SOtxaLJbEvtuNGzcCcOnSpcTZEwIDA4mNhUmToFQpGDjQHLmtVAl+/BE2bzZHbUVEREQcQeHWSbm5uVG+fHng7q0JFy5c4NKlS1gsFsqVK5dm17aH2w0bNgD/jtqWKlWKn37ypXx56NULQkLMUdvZs2HXLrO/9j5T8oqIiIikK4VbJ1a5cmUA9uzZc8fP7aO2xYsXv+8MBslRr149wBy5NQwj8WGyCxeq0rUrHD8OBQqYMx8cPGjOhODikmaXFxEREUkxPVDmxO43Y0JatyTYVatWDQ8PDy5dusQ778xi7NhPAAgP9yVnTujfH3r3Bm/vNL2siIiISKop3Doxe7i928itfY7bChUqpOl13d3dKVKkHEeP7mbMmBcTt9eo4c2qVZAjR5peTkRERCTNqC3BidnbEg4ePEhsbOxtn9vDbVqO3G7ZYi6Ve/SocWOLlSpVXmfDhkP8+edoBVsRERFxagq3Tqxo0aJkz56d+Ph4goODb/s8LUdud+2C1q2hdm1zqVyr9QPy5KnO1KkL2bnzS+rWLYunp2eqryMiIiKSnhRunZjFYrnrQ2Xh4eGcOXMGSF24PXAAnn0WAgNh6VKwWuGFF+DIkXZcuLCVl15qm+Jzi4iIiGQ0hVsnZw+3u3fvvmX7wYMHAShYsCB+fn7JPu/Ro9C1KwQEwA8/mFN4dewI+/fDzJnmFF8iIiIimY0eKHNydxu5TWlLwsmT8NFHMGMGJCSY29q2hWHD4MalRERERDIthVsnl1bh9tw5c4ncKVPA/mxaixYwfDjUqJF29YqIiIg4ksKtk7OH25MnTxIWFpbYgpDUcHvhAoweDRMmQHS0ua1xY/jwQ7ixVoOIiIhIlqGeWyeXM2dOihQpAty6mIO959a+RO9/Xb0K778PJUvCZ5+ZwbZOHQgKgjVrFGxFREQka1K4zQQeeugh4N+HymJjYzl69Chw+8hteLjZU1uihPnPiAioVg2WL4cNG6BJk4ytXURERCQjKdxmAv8Nt0eOHCEhIYHs2bNTsGBBwAy1I0aYofb9982R24AA+PFH2LrV7K+1WBz1DUREREQyhnpuM4H/hlt7v2358uW5ds3ChAkwZgxcvmzuX64cDB0KHTqY89aKiIiIPCgUbjMBe7jds2cPNpstsd82Pr4CJUrcGmqHDDFDrYuLo6oVERERcRyF20ygbNmyuLu7c+3aNbZtO8RXXy0BYMcOf0ChVkRERMRO4TYTcHNzo1y5iuzZs5NateoBVwDInv0Ykycr1IqIiIjYKdw6uWvXzDlq9+69cGPLFazWvNSs2ZjvvhtBqVIOLU9ERETEqSjcOqlr12DiRHOO2kuXAFoCs6hduy1r1szGy8vDwRWKiIiIOB89S+9krl2DUaPMKb0GDjSDbdmy8O23X3P9+jU2b56vYCsiIiJyF04ZbidOnEjx4sXx9PSkdu3abNmy5Z77//DDD5QvXx5PT08qV67M8uXLM6jStHOnUFumDHz7LezbB88/b1GoFREREbkPpwu38+fPp1+/fgwdOpTt27dTpUoVmjdvzvnz5++4/8aNG+nYsSM9evRgx44dtGnThjZt2tyyVK0zu1eo3b8fnn8eXNU8IiIiIpIkFsMwDEcXcbPatWtTs2ZNJkyYAIDNZsPf35/XX3+dAQMG3LZ/hw4duH79Or/88kvitocffpjAwEAmT5583+uFh4fj5+dHWFgYvr6+afdF7iMiwuyp/fRTe0+tGWrffx86dlSgFREREbFLTl5zqpHb2NhYtm3bRrNmzRK3Wa1WmjVrxqZNm+54zKZNm27ZH6B58+Z33T8mJobw8PBbXhkpIgI++QSKF4cBA/4dqZ092xyp7dJFwVZEREQkpZwq3F68eJGEhATy589/y/b8+fMTEhJyx2NCQkKStf/IkSPx8/NLfPn7+6dN8UnUps2/obZ0aYVaERERkbTkVOE2IwwcOJCwsLDE16lTpzL0+r16maH2m2/gwAGFWhEREZG05FSxKk+ePLi4uBAaGnrL9tDQUAoUKHDHYwoUKJCs/T08PPDwcNysA61bQ6tWCrQiIiIi6cGpRm7d3d2pXr06QUFBidtsNhtBQUHUqVPnjsfUqVPnlv0BVq1addf9Hc1qVbAVERERSS9OF7P69etHt27dqFGjBrVq1WLs2LFcv36d7t27A9C1a1cKFy7MyJEjAejTpw8NGzZkzJgxPPHEE8ybN4+tW7cyZcoUR34NEREREXEApwu3HTp04MKFCwwZMoSQkBACAwNZsWJF4kNjJ0+exGr9d8C5bt26zJ07l8GDB/Pee+9RpkwZlixZQkBAgKO+goiIiIg4iNPNc5vRHDXPrYiIiIgkTaad51ZEREREJDUUbkVEREQky1C4FREREZEsQ+FWRERERLIMhVsRERERyTIUbkVEREQky1C4FREREZEsQ+FWRERERLIMhVsRERERyTIUbkVEREQky1C4FREREZEsQ+FWRERERLIMhVsRERERyTIUbkVEREQky3B1dAGOZhgGAOHh4Q6uRERERETuxJ7T7LntXh74cHvt2jUA/P39HVyJiIiIiNzLtWvX8PPzu+c+FiMpETgLs9lsnD17luzZs2OxWBxdjkOFh4fj7+/PqVOn8PX1dXQ5Tk33Kul0r5JO9yppdJ+STvcq6XSvks4R98owDK5du0ahQoWwWu/dVfvAj9xarVaKFCni6DKciq+vr/7DTiLdq6TTvUo63auk0X1KOt2rpNO9SrqMvlf3G7G10wNlIiIiIpJlKNyKiIiISJahcCuJPDw8GDp0KB4eHo4uxenpXiWd7lXS6V4lje5T0uleJZ3uVdI5+7164B8oExEREZGsQyO3IiIiIpJlKNyKiIiISJahcCsiIiIiWYbCrYiIiIhkGQq3D5iJEydSvHhxPD09qV27Nlu2bLnrvosXL6ZGjRrkyJEDb29vAgMD+fbbbzOwWsdKzr262bx587BYLLRp0yZ9C3QiyblXs2bNwmKx3PLy9PTMwGodJ7m/pq5evUqvXr0oWLAgHh4elC1bluXLl2dQtY6VnHvVqFGj235NWSwWnnjiiQys2HGS++tq7NixlCtXjmzZsuHv70/fvn2Jjo7OoGodKzn3Ki4ujuHDh1OqVCk8PT2pUqUKK1asyMBqHWPdunW0atWKQoUKYbFYWLJkyX2PWbt2LdWqVcPDw4PSpUsza9asdK/zngx5YMybN89wd3c3ZsyYYezbt894+eWXjRw5chihoaF33P/33383Fi9ebOzfv984cuSIMXbsWMPFxcVYsWJFBlee8ZJ7r+yOHz9uFC5c2GjQoIHRunXrjCnWwZJ7r2bOnGn4+voa586dS3yFhIRkcNUZL7n3KSYmxqhRo4bRsmVLY/369cbx48eNtWvXGjt37szgyjNecu/VpUuXbvn1tHfvXsPFxcWYOXNmxhbuAMm9V999953h4eFhfPfdd8bx48eNlStXGgULFjT69u2bwZVnvOTeq3fffdcoVKiQsWzZMuPo0aPGpEmTDE9PT2P79u0ZXHnGWr58uTFo0CBj8eLFBmD8+OOP99z/2LFjhpeXl9GvXz9j//79xvjx4x2eFRRuHyC1atUyevXqlfg+ISHBKFSokDFy5Mgkn6Nq1arG4MGD06M8p5KSexUfH2/UrVvXmDZtmtGtW7cHJtwm917NnDnT8PPzy6DqnEdy79NXX31llCxZ0oiNjc2oEp1Gan+v+uKLL4zs2bMbERER6VWi00juverVq5fRpEmTW7b169fPqFevXrrW6QySe68KFixoTJgw4ZZt7dq1Mzp37pyudTqTpITbd99916hUqdIt2zp06GA0b948HSu7N7UlPCBiY2PZtm0bzZo1S9xmtVpp1qwZmzZtuu/xhmEQFBTEoUOHeOSRR9KzVIdL6b0aPnw4+fLlo0ePHhlRplNI6b2KiIigWLFi+Pv707p1a/bt25cR5TpMSu7T0qVLqVOnDr169SJ//vwEBAQwYsQIEhISMqpsh0jt71UA06dP57nnnsPb2zu9ynQKKblXdevWZdu2bYl/HX/s2DGWL19Oy5YtM6RmR0nJvYqJibmtZSpbtmysX78+XWvNbDZt2nTLfQVo3rx5kv97TQ+uDruyZKiLFy+SkJBA/vz5b9meP39+Dh48eNfjwsLCKFy4MDExMbi4uDBp0iQeffTR9C7XoVJyr9avX8/06dPZuXNnBlToPFJyr8qVK8eMGTN46KGHCAsL47PPPqNu3brs27ePIkWKZETZGS4l9+nYsWOsWbOGzp07s3z5co4cOcJrr71GXFwcQ4cOzYiyHSKlv1fZbdmyhb179zJ9+vT0KtFppORederUiYsXL1K/fn0MwyA+Pp6ePXvy3nvvZUTJDpOSe9W8eXM+//xzHnnkEUqVKkVQUBCLFy/O8n/ATK6QkJA73tfw8HCioqLIli1bhtekkVu5p+zZs7Nz507+/vtvPv74Y/r168fatWsdXZZTuXbtGl26dGHq1KnkyZPH0eU4vTp16tC1a1cCAwNp2LAhixcvJm/evHz99deOLs2p2Gw28uXLx5QpU6hevTodOnRg0KBBTJ482dGlObXp06dTuXJlatWq5ehSnNLatWsZMWIEkyZNYvv27SxevJhly5bx4YcfOro0pzNu3DjKlClD+fLlcXd3p3fv3nTv3h2rVdHJ2Wnk9gGRJ08eXFxcCA0NvWV7aGgoBQoUuOtxVquV0qVLAxAYGMiBAwcYOXIkjRo1Ss9yHSq59+ro0aOcOHGCVq1aJW6z2WwAuLq6cujQIUqVKpW+RTtISn9d3czNzY2qVaty5MiR9CjRKaTkPhUsWBA3NzdcXFwSt1WoUIGQkBBiY2Nxd3dP15odJTW/pq5fv868efMYPnx4epboNFJyr95//326dOnCSy+9BEDlypW5fv06r7zyCoMGDcqywS0l9ypv3rwsWbKE6OhoLl26RKFChRgwYAAlS5bMiJIzjQIFCtzxvvr6+jpk1BY0cvvAcHd3p3r16gQFBSVus9lsBAUFUadOnSSfx2azERMTkx4lOo3k3qvy5cuzZ88edu7cmfh66qmnaNy4MTt37sTf3z8jy89QafHrKiEhgT179lCwYMH0KtPhUnKf6tWrx5EjRxL/oAQQHBxMwYIFs2ywhdT9mvrhhx+IiYnh+eefT+8ynUJK7lVkZORtAdb+ByjDMNKvWAdLza8rT09PChcuTHx8PIsWLaJ169bpXW6m8v927hgU/j+O4/hHdN9v6YqEidydkkGx0BWdjAYjlus2gwwWdWVAcQtZZLLYXGSRxRnIoFt0d4PDcDrKrBhdvX6Tq//vN90/fT/18XyUQb7q/XlHPePuE4/H/7NXY4y5vLxsqC1+nLW3siFw2WxWnufp8PBQ5XJZCwsLamtrq1/DlEwmlU6n689nMhnlcjlVKhWVy2Xt7OyopaVFBwcHto4QmEZ39bffdFtCo7va2NjQxcWFKpWK7u7uND8/L9/3dX9/b+sIgWh0T6+vrwqHw1paWtLT05POz8/V1dWlzc1NW0cIzP/9/RsfH9fc3FzQ41rV6K7W1tYUDod1dHSk5+dn5XI5xWIxzc7O2jpCYBrdVT6f1+npqSqVim5ubjQ1NaVIJKL393dLJwjG5+enCoWCCoWCjDHa3d1VoVDQy8uLJCmdTiuZTNaf/74KbGVlRQ8PD9rf3+cqMARrb29Pvb29CoVCGh0dVT6fr38tkUgolUrVP19dXVV/f79831d7e7vi8biy2ayFqe1oZFd/+01xKzW2q+Xl5fqz3d3dmp6edv7eyG+N/kzd3t5qbGxMnucpGo1qa2tLtVot4KntaHRXj4+PMsYol8sFPKl9jezq6+tL6+vrisVi8n1fPT09WlxcdD7YvjWyq+vraw0ODsrzPHV0dCiZTOrt7c3C1MG6urqSMeafj+/dpFIpJRKJf75neHhYoVBI0WjU+h3TTZLD/4cAAADAr8JrbgEAAOAM4hYAAADOIG4BAADgDOIWAAAAziBuAQAA4AziFgAAAM4gbgEAAOAM4hYAAADOIG4BAADgDOIWABxXq9VsjwAAgSFuAcAh1WrVNDU1mePjYzMxMWE8zzNnZ2e2xwKAwLTYHgAA8HNKpZIxxpjt7W2TyWRMJBIxnZ2dlqcCgOAQtwDgkGKxaFpbW83JyYnp6+uzPQ4ABI6XJQCAQ0qlkpmZmSFsAfxaxC0AOKRYLJrJyUnbYwCANcQtADji4+PDVKtVMzIyYnsUALCGuAUAR5RKJdPc3GyGhoZsjwIA1hC3AOCIUqlkBgYGjO/7tkcBAGuaJMn2EAAAAMBP4C+3AAAAcAZxCwAAAGcQtwAAAHAGcQsAAABnELcAAABwBnELAAAAZxC3AAAAcAZxCwAAAGcQtwAAAHAGcQsAAABnELcAAABwBnELAAAAZ/wBBDY1lTadMSYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot homogenised potential\n",
    "fig, ax = plt.subplots(1, 1, figsize=(8, 6))\n",
    "\n",
    "ax.plot(r_total_mesh, phi_n(r=r_total_mesh), \"b\", label=r\"$\\phi^-$\")\n",
    "ax.plot(r_total_mesh, phi_p(r=r_total_mesh), \"r\", label=r\"$\\phi^+$\")\n",
    "for i in range(len(tt)):\n",
    "    ax.plot(\n",
    "        r_mesh_pos[i, :],\n",
    "        phi_p(r=r_mesh_pos[i, :]),\n",
    "        \"k\",\n",
    "        label=r\"$\\phi$\" if i == 0 else \"\",\n",
    "    )\n",
    "for i in range(len(tt) - 1):\n",
    "    ax.plot(r_mesh_neg[i, :], phi_n(r=r_mesh_neg[i, :]), \"k\")\n",
    "    ax.plot(r_mesh_am1[i, :], phi_am1(r_mesh_am1[i, :], tt[i]), \"k\")\n",
    "    ax.plot(r_mesh_am2[i, :], phi_am2(r_mesh_am2[i, :], tt[i]), \"k\")\n",
    "ax.set_xlabel(r\"$r$\")\n",
    "ax.set_ylabel(r\"$\\phi$\")\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "thrown-proposition",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "The relevant papers for this notebook are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "vietnamese-brisbane",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, 2019. doi:10.1007/s12532-018-0139-4.\n",
      "[2] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.\n",
      "[3] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "pybamm.print_citations()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.6"
  },
  "vscode": {
   "interpreter": {
    "hash": "9ff3d0c7e37de5f5aa47f4f719e4c84fc6cba7b39c571a05173422444e82fa58"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
